BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EsTimeAlg.cxx
Go to the documentation of this file.
1#include <algorithm>
2#include <iostream>
3#include <vector>
4
5#include "Emc_helix.h"
6#include "EsTimeAlg.h"
7#include "EstParameter.h"
8#include "Toffz_helix.h"
9#include "TrackUtil/Helix.h"
10
11#include "AIDA/IHistogramFactory.h"
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/IDataManagerSvc.h"
14#include "GaudiKernel/IDataProviderSvc.h"
15#include "GaudiKernel/IHistogramSvc.h"
16#include "GaudiKernel/IMessageSvc.h"
17#include "GaudiKernel/INTupleSvc.h"
18#include "GaudiKernel/ISvcLocator.h"
19#include "GaudiKernel/MsgStream.h"
20#include "GaudiKernel/PropertyMgr.h"
21#include "GaudiKernel/SmartDataPtr.h"
22#include "GaudiKernel/StatusCode.h"
23
24#include "CLHEP/Vector/ThreeVector.h"
25#include "DetVerSvc/IDetVerSvc.h"
26#include "EmcRecEventModel/RecEmcShower.h"
27#include "EstTofCaliSvc/IEstTofCaliSvc.h"
28#include "EventModel/EventHeader.h"
29#include "GaudiKernel/IPartPropSvc.h"
30#include "McTruth/McParticle.h"
31#include "McTruth/TofMcHit.h"
32#include "MdcRawEvent/MdcDigi.h"
33#include "MdcRecEvent/RecMdcDedx.h"
34#include "MdcRecEvent/RecMdcHit.h"
35#include "MdcRecEvent/RecMdcTrack.h"
36#include "RawDataProviderSvc/IRawDataProviderSvc.h"
37#include "RawDataProviderSvc/TofData.h"
38#include "RawEvent/RawDataUtil.h"
39#include "ReconEvent/ReconEvent.h"
40#include "TofRawEvent/TofDigi.h"
41#include "TrigEvent/TrigData.h"
42
43#include "MdcGeomSvc/IMdcGeomSvc.h"
44#include "MdcGeomSvc/MdcGeoLayer.h"
45#include "MdcGeomSvc/MdcGeoWire.h"
46// #include "MdcCalibFunSvc/MdcCalibFunSvc.h"
47#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
48
49#include "EvTimeEvent/RecEsTime.h"
50#include "Identifier/Identifier.h"
51#include "Identifier/MdcID.h"
52#include "Identifier/TofID.h"
53// #include "MdcFastTrkAlg/ntupleItem.h"
54
55#include "CLHEP/Geometry/Point3D.h"
56#include "CLHEP/Vector/ThreeVector.h"
57#ifndef ENABLE_BACKWARDS_COMPATIBILITY
58typedef HepGeom::Point3D<double> HepPoint3D;
59#endif
60
61using CLHEP::Hep3Vector;
62using CLHEP::HepMatrix;
63using CLHEP::HepSymMatrix;
64using CLHEP::HepVector;
65typedef std::vector<double> Vdouble;
66
67using namespace std;
68using namespace Event;
69
70// Constants
71
72const double VLIGHT = 29.98;
73const double ELMAS2 = 0.511E-3 * 0.511E-3;
74const double MUMAS2 = 105.658E-3 * 105.658E-3;
75const double PIMAS2 = 139.569E-3 * 139.569E-3;
76const double PROTONMAS2 = 938.272E-3 * 938.272E-3;
77const double RCTOF2 = 81. * 81.;
78const double RCEMC2 = 94.2 * 94.2;
79const double TDC2NSEC = 0.5, NSEC2TDC = 2.0;
80const double PI = 3.141593, PIBY44 = 3.141593 / 44.;
81const double VELPROP = 33.33;
82DECLARE_COMPONENT( EsTimeAlg )
83EsTimeAlg::EsTimeAlg( const std::string& name, ISvcLocator* pSvcLocator )
84 : Algorithm( name, pSvcLocator ) {
85 for ( int i = 0; i < 5; i++ ) m_pass[i] = 0;
86 m_flag = 1;
87 m_nbunch = 3;
88 m_ntupleflag = 1;
89 m_beforrec = 1;
90 m_optCosmic = 0;
91 m_cosmicScheme = 0;
92 m_evtNo = 0;
93 m_ebeam = 1.85;
94 m_userawtime_opt = 0;
95 m_debug = 0;
96 declareProperty( "MdcMethod", m_flag );
97 declareProperty( "Nbunch", m_nbunch );
98 declareProperty( "BunchtimeMC",
99 m_bunchtime_MC = 8.0 ); // assign bunch interval for MC; for data it's
100 // always obtained from calibration constants.
101 declareProperty( "NtupleFlag", m_ntupleflag );
102 declareProperty( "beforrec", m_beforrec );
103 declareProperty( "Cosmic", m_optCosmic );
104 declareProperty( "CosmScheme", m_cosmicScheme );
105 declareProperty( "EventNo", m_evtNo );
106 declareProperty( "Ebeam", m_ebeam );
107 declareProperty( "UseRawTime", m_userawtime_opt );
108 declareProperty( "RawOffset_b", toffset_rawtime = 0.0 );
109 declareProperty( "RawOffset_e", toffset_rawtime_e = 0.0 );
110 declareProperty( "Offset_dt_b", offset_dt = 0.0 );
111 declareProperty( "Offset_dt_e", offset_dt_end = 0.0 );
112 declareProperty( "debug", m_debug );
113 declareProperty( "UseXT", m_useXT = true );
114 declareProperty( "UseT0", m_useT0cal = true );
115 declareProperty( "UseSw", m_useSw = false );
116 declareProperty( "MdcOpt", m_mdcopt = false );
117 declareProperty( "TofOpt", m_TofOpt = false );
118 declareProperty( "TofOptCut", m_TofOpt_Cut = 12. );
119 declareProperty( "UseTimeFactor", m_useTimeFactor = true );
120}
121
122StatusCode EsTimeAlg::initialize() {
123
124 MsgStream log( msgSvc(), name() );
125 log << MSG::INFO << "in initialize()" << endmsg;
126 if ( m_ntupleflag )
127 {
128
129 NTuplePtr nt2( ntupleSvc(), "FILE105/h2" );
130
131 if ( nt2 ) m_tuple2 = nt2;
132 else
133 {
134 m_tuple2 = ntupleSvc()->book( "FILE105/h2", CLID_ColumnWiseTuple, "Event Start Time" );
135
136 if ( m_tuple2 )
137 {
138
139 m_tuple2->addItem( "eventNo", g_eventNo );
140 m_tuple2->addItem( "runNo", g_runNo );
141 // #ifndef McTruth
142 m_tuple2->addItem( "NtrackMC", g_ntrkMC, 0, 5000 );
143 m_tuple2->addItem( "MCtheta0", g_ntrkMC, g_theta0MC );
144 m_tuple2->addItem( "MCphi0", g_ntrkMC, g_phi0MC );
145 m_tuple2->addItem( "pxMC", g_ntrkMC, g_pxMC );
146 m_tuple2->addItem( "pyMC", g_ntrkMC, g_pyMC );
147 m_tuple2->addItem( "pzMC", g_ntrkMC, g_pzMC );
148 m_tuple2->addItem( "ptMC", g_ntrkMC, g_ptMC );
149 m_tuple2->addItem( "mct0", g_mcTestime );
150 // #endif
151 m_tuple2->addItem( "Ntrack", g_ntrk, 0, 5000 );
152 m_tuple2->addItem( "ttof", g_ntrk, g_ttof );
153 m_tuple2->addItem( "velocity", g_ntrk, g_vel );
154 m_tuple2->addItem( "abmom", g_ntrk, g_abmom );
155 m_tuple2->addItem( "pid", g_ntrk, g_pid );
156 m_tuple2->addItem( "nmatchBarrel", g_nmatchbarrel );
157 m_tuple2->addItem( "nmatchBarrel_1", g_nmatchbarrel_1 );
158 m_tuple2->addItem( "nmatchBarrel_2", g_nmatchbarrel_2 );
159 m_tuple2->addItem( "nmatchend", g_nmatchend );
160 m_tuple2->addItem( "nmatch", g_nmatch_tot );
161 m_tuple2->addItem( "t0forward", g_ntrk, g_t0for );
162 m_tuple2->addItem( "t0backward", g_ntrk, g_t0back );
163 m_tuple2->addItem( "meant0", g_meant0 );
164 m_tuple2->addItem( "nmatchmdc", g_nmatchmdc );
165 m_tuple2->addItem( "ndriftt", g_ndriftt );
166 m_tuple2->addItem( "MdcEsTime", g_EstimeMdc );
167 m_tuple2->addItem( "Mdct0mean", g_t0mean );
168 m_tuple2->addItem( "Mdct0try", g_t0 );
169 m_tuple2->addItem( "Mdct0sq", g_T0 );
170 m_tuple2->addItem( "trigtiming", g_trigtiming );
171 m_tuple2->addItem( "meantdc", g_meantdc );
172 // m_tuple2->addItem ( "meantev" , g_meantev);
173 m_tuple2->addItem( "ntofup", g_ntofup, 0, 500 );
174 m_tuple2->addItem( "ntofdown", g_ntofdown, 0, 500 );
175 m_tuple2->addIndexedItem( "meantevup", g_ntofup, g_meantevup );
176 m_tuple2->addIndexedItem( "meantevdown", g_ntofdown, g_meantevdown );
177 m_tuple2->addItem( "ntofup1", g_ntofup1 );
178 m_tuple2->addItem( "ntofdown1", g_ntofdown1 );
179 m_tuple2->addItem( "Testime_fzisan", g_Testime_fzisan );
180 m_tuple2->addItem( "Testime", g_Testime );
181 m_tuple2->addItem( "T0barrelTof", g_t0barrelTof );
182 m_tuple2->addItem( "difftofb", g_difftof_b );
183 m_tuple2->addItem( "difftofe", g_difftof_e );
184 m_tuple2->addItem( "EstFlag", m_estFlag );
185 m_tuple2->addItem( "EstTime", m_estTime );
186 }
187 else
188 { // did not manage to book the N tuple....
189 log << MSG::ERROR << "Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
190 }
191 }
192 NTuplePtr nt9( ntupleSvc(), "FILE105/h9" );
193
194 if ( nt9 ) m_tuple9 = nt9;
195 else
196 {
197 m_tuple9 = ntupleSvc()->book( "FILE105/h9", CLID_ColumnWiseTuple, "Event Start time" );
198
199 if ( m_tuple9 )
200 {
201 m_tuple9->addItem( "Nmatch", g_nmatch, 0, 500 );
202 m_tuple9->addIndexedItem( "meantev", g_nmatch, g_meantev );
203 }
204 else { log << MSG::ERROR << "Cannot book N-tuple:" << long( m_tuple9 ) << endmsg; }
205 }
206 NTuplePtr nt3( ntupleSvc(), "FILE105/calibconst" );
207
208 if ( nt3 ) m_tuple3 = nt3;
209 else
210 {
211 m_tuple3 =
212 ntupleSvc()->book( "FILE105/calibconst", CLID_ColumnWiseTuple, "Event Start time" );
213
214 if ( m_tuple3 )
215 {
216 m_tuple3->addItem( "t0offsetb", g_t0offset_b );
217 m_tuple3->addItem( "t0offsete", g_t0offset_e );
218 m_tuple3->addItem( "bunchtime", g_bunchtime );
219 }
220 else { log << MSG::ERROR << "Cannot book N-tuple:" << long( m_tuple3 ) << endmsg; }
221 }
222 }
223
224 // Get the Particle Properties Service
225 IPartPropSvc* p_PartPropSvc;
226 static const bool CREATEIFNOTTHERE( true );
227 StatusCode PartPropStatus = service( "PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE );
228 if ( !PartPropStatus.isSuccess() || 0 == p_PartPropSvc )
229 {
230 log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
231 return PartPropStatus;
232 }
233 m_particleTable = p_PartPropSvc->PDT();
234 // IRawDataProviderSvc* m_rawDataProviderSvc;
235 StatusCode RawDataStatus =
236 service( "RawDataProviderSvc", m_rawDataProviderSvc, CREATEIFNOTTHERE );
237 if ( !RawDataStatus.isSuccess() )
238 {
239 log << MSG::ERROR << "Could not load RawDataProviderSvc!" << m_rawDataProviderSvc
240 << endmsg;
241 return RawDataStatus;
242 }
243
244 IDetVerSvc* detVerSvc;
245 StatusCode sc_det = service( "DetVerSvc", detVerSvc );
246 if ( sc_det.isFailure() )
247 {
248 log << MSG::ERROR << "can't retrieve DetVerSvc instance" << endmsg;
249 return sc_det;
250 }
251 m_phase = detVerSvc->phase();
252
253 StatusCode sc = service( "CalibDataSvc", m_pCalibDataSvc, true );
254 if ( !sc.isSuccess() )
255 {
256 log << MSG::ERROR << "Could not get IDataProviderSvc interface of CalibXmlCnvSvc"
257 << endmsg;
258 return sc;
259 }
260 else
261 { log << MSG::DEBUG << "Retrieved IDataProviderSvc interface of CalibXmlCnvSvc" << endmsg; }
262 /*maqm comment
263 sc = service("CalibRootCnvSvc", m_pRootSvc, true);
264 if ( !sc.isSuccess() ) {
265 log << MSG::ERROR
266 << "Could not get ICalibRootSvc interface of CalibRootCnvSvc"
267 << endmsg;
268 return sc;
269 }*/
270
271 // Get properties from the JobOptionsSvc
272
273 // sc = setProperties();
274
275 // Get TOF Calibtration Service
276 // IEstTofCaliSvc* tofCaliSvc;
277 StatusCode scc = service( "EstTofCaliSvc", tofCaliSvc );
278 if ( scc == StatusCode::SUCCESS )
279 {
280 // log<<MSG::INFO<<"begin dump Calibration Constants"<<endmsg;
281 // tofCaliSvc->Dump();
282 log << MSG::INFO << " Get EstTof Calibration Service Sucessfully!! " << endmsg;
283 }
284 else
285 {
286 log << MSG::ERROR << " Get EstTof Calibration Service Failed !! " << endmsg;
287 return scc;
288 }
289
290 if ( m_useXT || m_useT0cal )
291 {
292 StatusCode scc = Gaudi::svcLocator()->service( "MdcCalibFunSvc", m_mdcCalibFunSvc );
293 if ( scc.isFailure() ) { log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endmsg; }
294 }
295
296 // Initailize MdcUtilitySvc
297
298 sc = service( "MdcUtilitySvc", m_mdcUtilitySvc );
299 if ( sc.isFailure() )
300 {
301 log << MSG::FATAL << "Could not load MdcUtilitySvc!" << endmsg;
302 return StatusCode::FAILURE;
303 }
304
305 return StatusCode::SUCCESS;
306}
307
308StatusCode EsTimeAlg::execute() {
309
310 MsgStream log( msgSvc(), name() );
311 log << MSG::INFO << " tof " << endmsg;
312
313 // Constants
314 static const EstParameter Estparam;
315 double offset = 0, t_quality = 0, tOffset_b = 0, tOffset_e = 0;
316 int idtof, tofid_helix[30] = { -9 }, idmatch[3][88] = { 0 }, idmatch_emc[3][88] = { 0 },
317 idt[15] = { 0 }, particleId[30] = { 0 }, tofid_emc[2] = { 0 }, module[20] = { 0 };
318 int idetf, etfid_helix[30] = { -9 }, idetfmatch[3][36] = { -9 },
319 idmatch_etf_emc[3][36] = { 0 }, etfid_emc[2] = { 0 };
320 int ntot = 0, in = -1, out = -1, emcflag1 = 0, emcflag2 = 0, tof_flag = 0;
321 double bunchtime = m_bunchtime_MC;
322 double meant[15] = { 0. }, adc[15] = { 0. }, momentum[15] = { 0. }, r_endtof[15] = { 0. };
323 double ttof[30] = { 0. }, helztof[30] = { 0.0 }, mcztof = 0.0, forevtime = 0.0,
324 backevtime = 0.0, meantev[500] = { 0. }, meantevup[500] = { 0.0 },
325 meantevdown[500] = { 0.0 };
326 double t0forward = 0, t0backward = 0, t0forward_trk = 0, t0backward_trk = 0;
327 double t0forward_add = 0, t0backward_add = 0, t_Est = -999;
328 double thetaemc_rec[20] = { 0. }, phiemc_rec[20] = { 0. }, energy_rec[20] = { 0. },
329 xemc_rec[20] = { 0. }, yemc_rec[20] = { 0. }, zemc_rec[20] = { 0. };
330 double r_endetf[30] = { 0. }, tetf[30] = { 0. }, helzetf[30] = { 0. },
331 helpathetf[36] = { 0. }, abmom2etf[36] = { 0. };
332
333 int nmatch1 = 0, nmatch2 = 0, nmatch_barrel = 0, nmatch_end = 0, nmatch_mdc = 0,
334 nmatch_barrel_1 = 0, nmatch_barrel_2 = 0, nmatch = 0, ntofup = 0, ntofdown = 0;
335 double sum_EstimeMdc = 0, sum_EstimeMdcMC = 0;
336 int nbunch = 0, tEstFlag = 0, runNo = 0;
337 double helpath[88] = { 0. }, helz[88] = { 0. }, abmom2[88] = { 0. };
338 double mcTestime = 0, trigtiming = 0;
339 double mean_tdc_btof[2][88] = { 0. }, mean_tdc_etof[3][48] = { 0. },
340 mean_tdc_etf[3][36][12] = { 0. };
341 int optCosmic = m_optCosmic;
342 int m_userawtime = m_userawtime_opt;
343 double Testime_fzisan = -999.; // yzhang add
344 int TestimeFlag_fzisan = -999; // yzhang add
345 double TestimeQuality_fzisan = -999.; // yzhang add
346 double Tof_t0Arr[500] = { -999. };
347
348 bool useEtofScin = ( m_phase < 3 );
349 bool useEtofMRPC = ( m_phase > 2 );
350
351 // Part 1: Get the event header, print out event and run number
352 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
353 if ( !eventHeader )
354 {
355 log << MSG::FATAL << "Could not find Event Header" << endmsg;
356 return StatusCode::FAILURE;
357 }
358 int eventNo = eventHeader->eventNumber();
359 runNo = eventHeader->runNumber();
360 log << MSG::INFO << "EsTimeAlg: retrieved event: " << eventNo << " run: " << runNo
361 << endmsg;
362
363 if ( m_ntupleflag && m_tuple2 )
364 {
365 g_runNo = runNo;
366 g_eventNo = eventNo;
367 }
368 if ( runNo < 0 )
369 {
370 offset_dt = 0;
371 offset_dt_end = 0;
372 toffset_rawtime = 0;
373 toffset_rawtime_e = 0;
374 if ( useEtofMRPC ) { toffset_rawtime_e = -20.0; }
375 }
376 if ( runNo > 0 && useEtofMRPC ) { toffset_rawtime_e = 1.7; }
377
378 if ( runNo > 0 && runNo < 5336 ) offset_dt = 5.8;
379 if ( runNo > 5462 && runNo < 5466 )
380 {
381 offset_dt = 1.6;
382 offset_dt_end = 0.4;
383 }
384 if ( runNo > 5538 && runNo < 5920 )
385 {
386 offset_dt = -0.2;
387 offset_dt_end = -1.2;
388 }
389 if ( runNo > 5924 && runNo < 6322 )
390 {
391 offset_dt = -0.2;
392 offset_dt_end = -0.1;
393 }
394 if ( runNo > 6400 && runNo < 6423 )
395 {
396 offset_dt = -2.4;
397 offset_dt_end = -1.9;
398 }
399 if ( runNo > 6514 && runNo < 6668 )
400 {
401 offset_dt = 0;
402 offset_dt_end = 3.4;
403 }
404 if ( runNo > 6667 && runNo < 6715 )
405 {
406 offset_dt = 4;
407 offset_dt_end = 7.2;
408 }
409 if ( runNo > 6715 && runNo < 6720 )
410 {
411 offset_dt = 8;
412 offset_dt_end = 7.5;
413 }
414 if ( runNo > 6879 && runNo < 6909 )
415 {
416 offset_dt = 0.2;
417 offset_dt_end = -0.1;
418 }
419 if ( m_ntupleflag && m_tuple3 )
420 {
421 g_bunchtime = 8;
422 g_t0offset_b = 2.0;
423 g_t0offset_e = 2.0;
424 }
425
426 m_pass[0]++;
427 if ( m_evtNo == 1 && m_pass[0] % 1000 == 0 )
428 { cout << "------------------- Events-----------------: " << m_pass[0] << endl; }
429 if ( m_debug == 4 ) cout << "m_userawtime: " << m_userawtime << endl;
430 if ( m_debug == 4 ) cout << "EstTofCalibSvc est flag: " << tofCaliSvc->ValidInfo() << endl;
431
432 if ( tofCaliSvc->chooseConstants( runNo, eventNo ) == StatusCode::FAILURE )
433 {
434 log << MSG::ERROR << "EstTof Calibration Const is NOT correct! " << endmsg;
435 return StatusCode::FAILURE;
436 }
437 if ( tofCaliSvc->ValidInfo() == 0 && m_userawtime_opt == 0 )
438 {
439 log << MSG::ERROR << "EstTof Calibration Const is Invalid! " << endmsg;
440 return StatusCode::FAILURE;
441 }
442 if ( tofCaliSvc->ValidInfo() == 0 || m_userawtime_opt == 1 ) m_userawtime = 1;
443 else m_userawtime = 0;
444
445 SmartDataPtr<RecEsTimeCol> aRecestimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
446 if ( !aRecestimeCol || aRecestimeCol->size() == 0 )
447 {
448 if ( m_debug == 4 )
449 log << MSG::INFO << "Could not find RecEsTimeCol from fzsian" << endmsg;
450 }
451 else
452 {
453 RecEsTimeCol::iterator it_evt = aRecestimeCol->begin();
454 for ( ; it_evt != aRecestimeCol->end(); it_evt++ )
455 {
456 Testime_fzisan = ( *it_evt )->getTest(); // yzhang add
457 TestimeFlag_fzisan = ( *it_evt )->getStat(); // yzhang add
458 TestimeQuality_fzisan = ( *it_evt )->getQuality(); // yzhang add
459
460 log << MSG::INFO << "fzisan : Test = " << ( *it_evt )->getTest()
461 << ", Status = " << ( *it_evt )->getStat() << endmsg;
462
463 if ( m_ntupleflag && m_tuple2 ) g_Testime_fzisan = Testime_fzisan;
464 }
465 }
466
467 static std::string fullPath = "/Calib/EsTimeCal";
468 SmartDataPtr<CalibData::EsTimeCalibData> TEst( m_pCalibDataSvc, fullPath );
469 if ( !TEst ) { cout << "ERROR EsTimeCalibData" << endl; }
470 else
471 {
472 int no = TEst->getTestCalibConstNo();
473 // std::cout<<"no==========="<<no<<std::endl;
474 // for(int i=0;i<no;i++){
475 // std::cout<<"i= "<<i<<" value="<< TEst->getTEstCalibConst(i)<<std::endl;
476 // }
477
478 unsigned int inumber = 0;
479 unsigned int calibNo = TEst->getSize();
480 if ( calibNo > 1 )
481 {
482 for ( unsigned int i = 0; i < calibNo; i++, inumber++ )
483 {
484 if ( ( TEst->getRunTo( i ) != -1 ) && ( TEst->getRunTo( i ) < TEst->getRunFrom( i ) ) )
485 {
486 log << MSG::ERROR << "EsTimeCal -- The " << inumber
487 << "th calibration constatns is ABNORMAL! Run From is LARGER than RUN To!"
488 << endmsg;
489 return StatusCode::FAILURE;
490 }
491 if ( ( TEst->getRunFrom( i ) == TEst->getRunTo( i ) ) &&
492 ( TEst->getEventFrom( i ) != -1 ) && ( TEst->getEventTo( i ) != -1 ) &&
493 ( TEst->getEventFrom( i ) > TEst->getEventTo( i ) ) )
494 {
495 log << MSG::ERROR << "EsTimeCal -- The " << inumber
496 << "th calibration constatns is ABNORMAL! Event From is LARGER than Event To!"
497 << endmsg;
498 return StatusCode::FAILURE;
499 }
500 }
501
502 inumber = 0;
503 bool filled = false;
504 for ( unsigned int i = 0; i < calibNo; i++, inumber++ )
505 {
506 int runFrom = TEst->getRunFrom( i );
507 int runTo = TEst->getRunTo( i );
508 int eventFrom = TEst->getEventFrom( i );
509 int eventTo = TEst->getEventTo( i );
510 if ( ( runNo == runFrom ) && ( ( eventFrom == -1 ) || ( eventNo >= eventFrom ) ) )
511 {
512 if ( ( runNo < runTo ) ||
513 ( ( runNo == runTo ) && ( ( eventTo == -1 ) || ( eventNo <= eventTo ) ) ) )
514 {
515 filled = true;
516 break;
517 }
518 }
519 if ( runNo > runFrom )
520 {
521 if ( ( runNo < runTo ) ||
522 ( ( runNo == runTo ) && ( ( eventTo == -1 ) || ( eventNo <= eventTo ) ) ) )
523 {
524 filled = true;
525 break;
526 }
527 }
528 }
529 if ( !filled )
530 {
531 log << MSG::ERROR << "EsTimeCal -- For run number:" << runNo
532 << ", NO suitable calibration constant is found!" << endmsg;
533 return StatusCode::FAILURE;
534 }
535 }
536
537 log << MSG::INFO << "offset barrel t0=" << TEst->getToffsetb( inumber )
538 << ", offset endcap t0=" << TEst->getToffsete( inumber )
539 << ", bunch time =" << TEst->getBunchTime( inumber ) << endmsg;
540 tOffset_b = TEst->getToffsetb( inumber );
541 tOffset_e = TEst->getToffsete( inumber );
542 bunchtime = TEst->getBunchTime( inumber );
543 }
544
545 if ( m_userawtime )
546 {
547 tOffset_b = 0;
548 tOffset_e = 0;
549 }
550 else
551 {
552 toffset_rawtime = 0;
553 toffset_rawtime_e = 0;
554 offset_dt = 0;
555 offset_dt_end = 0;
556 }
557
558 if ( runNo > 0 && m_useTimeFactor )
559 { bunchtime = RawDataUtil::TofTime( 8 * 1024 ) * bunchtime / 24.0; }
560
561 if ( runNo < 0 ) bunchtime = m_bunchtime_MC;
562
563 // Part 2: Retrieve MC truth
564 int digiId;
565 if ( runNo < 0 )
566 {
567 SmartDataPtr<McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
568 if ( !mcParticleCol ) { log << MSG::INFO << "Could not find McParticle" << endmsg; }
569 else
570 {
571 McParticleCol::iterator iter_mc = mcParticleCol->begin();
572 digiId = 0;
573 ntrkMC = 0;
574 int charge = 0;
575
576 for ( ; iter_mc != mcParticleCol->end(); iter_mc++, digiId++ )
577 {
578 int statusFlags = ( *iter_mc )->statusFlags();
579 int pid = ( *iter_mc )->particleProperty();
580 log << MSG::INFO << " MC ParticleId = " << pid << " statusFlags = " << statusFlags
581 << " PrimaryParticle = " << ( *iter_mc )->primaryParticle() << endmsg;
582 if ( m_ntupleflag && m_tuple2 )
583 {
584 g_theta0MC[ntrkMC] = ( *iter_mc )->initialFourMomentum().theta();
585 g_phi0MC[ntrkMC] = ( *iter_mc )->initialFourMomentum().phi();
586 g_pxMC[ntrkMC] = ( *iter_mc )->initialFourMomentum().px() / 1000;
587 g_pyMC[ntrkMC] = ( *iter_mc )->initialFourMomentum().py() / 1000;
588 g_pzMC[ntrkMC] = ( *iter_mc )->initialFourMomentum().pz() / 1000;
589 g_ptMC[ntrkMC] = sqrt( ( ( *iter_mc )->initialFourMomentum().px() ) *
590 ( ( *iter_mc )->initialFourMomentum().px() ) +
591 ( ( *iter_mc )->initialFourMomentum().py() ) *
592 ( ( *iter_mc )->initialFourMomentum().py() ) ) /
593 1000;
594 }
595 if ( pid > 0 )
596 {
597 if ( m_particleTable->particle( pid ) )
598 charge = m_particleTable->particle( pid )->charge();
599 }
600 else if ( pid < 0 )
601 {
602 if ( m_particleTable->particle( -pid ) )
603 {
604 charge = m_particleTable->particle( -pid )->charge();
605 charge *= -1;
606 }
607 }
608 else { log << MSG::WARNING << "wrong particle id, please check data" << endmsg; }
609 log << MSG::DEBUG << "MC ParticleId = " << pid << " charge = " << charge << endmsg;
610 if ( charge != 0 && abs( cos( ( *iter_mc )->initialFourMomentum().theta() ) ) < 0.93 )
611 { ntrkMC++; }
612 if ( ( ( *iter_mc )->primaryParticle() ) && m_ntupleflag && m_tuple2 )
613 {
614 g_mcTestime = ( *iter_mc )->initialPosition().t();
615 mcTestime = ( *iter_mc )->initialPosition().t();
616 }
617 }
618 if ( m_ntupleflag && m_tuple2 ) g_ntrkMC = ntrkMC;
619 }
620 } // close if(runno<0)
621 if ( m_debug ) cout << "bunchtime: " << bunchtime << endl;
622
623 SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc(), "/Event/Recon/RecMdcTrackCol" );
624 if ( !newtrkCol || newtrkCol->size() == 0 )
625 { log << MSG::INFO << "Could not find RecMdcTrackCol" << endmsg; }
626 else
627 {
628 log << MSG::INFO << "Begin to check RecMdcTrackCol" << endmsg;
629 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
630 for ( ; iter_trk != newtrkCol->end(); iter_trk++ )
631 {
632 log << MSG::DEBUG << "retrieved MDC track:"
633 << " Track Id: " << ( *iter_trk )->trackId()
634 << " Phi0: " << ( *iter_trk )->helix( 0 ) << " kappa: " << ( *iter_trk )->helix( 2 )
635 << " Tanl: " << ( *iter_trk )->helix( 4 )
636 << " Phi terminal: " << ( *iter_trk )->getFiTerm() << endmsg
637 << "Number of hits: " << ( *iter_trk )->getNhits() << " Number of stereo hits "
638 << ( *iter_trk )->nster() << endmsg;
639 double kappa = ( *iter_trk )->helix( 2 );
640 double tanl = ( *iter_trk )->helix( 4 );
641 if ( ( *iter_trk )->helix( 3 ) > 50.0 ) continue;
642 ntot++;
643 if ( ntot > 14 ) break;
644 momentum[ntot] = 1. / fabs( kappa ) * sqrt( 1. + tanl * tanl );
645 }
646 }
647
648 // get EmcRec results
649 int emctrk = 0;
650 SmartDataPtr<RecEmcShowerCol> aShowerCol( eventSvc(), "/Event/Recon/RecEmcShowerCol" );
651 if ( !aShowerCol || aShowerCol->size() == 0 )
652 { log << MSG::WARNING << "Could not find RecEmcShowerCol" << endmsg; }
653 else
654 {
655 RecEmcShowerCol::iterator iShowerCol = aShowerCol->begin();
656 for ( ; iShowerCol != aShowerCol->end(); iShowerCol++, emctrk++ )
657 {
658 if ( emctrk > 19 ) break;
659 phiemc_rec[emctrk] = ( *iShowerCol )->position().phi();
660 thetaemc_rec[emctrk] = ( *iShowerCol )->position().theta();
661 energy_rec[emctrk] = ( *iShowerCol )->energy();
662 xemc_rec[emctrk] = ( *iShowerCol )->x();
663 yemc_rec[emctrk] = ( *iShowerCol )->y();
664 zemc_rec[emctrk] = ( *iShowerCol )->z();
665 module[emctrk] = ( *iShowerCol )->module();
666 }
667 }
668 for ( int i = 0; i < 2; i++ )
669 {
670 double fi_endtof = atan2( yemc_rec[i], xemc_rec[i] ); // atna2 from -pi to pi
671 if ( fi_endtof < 0 ) { fi_endtof = 2 * 3.141593 + fi_endtof; }
672 if ( module[i] == 1 )
673 {
674 int Tofid = (int)( fi_endtof / ( 3.141593 / 44 ) );
675 if ( Tofid > 87 ) Tofid = Tofid - 88;
676 tofid_emc[i] = Tofid;
677 idmatch_emc[1][Tofid] = 1;
678 }
679 else
680 {
681 if ( useEtofScin )
682 {
683 int Tofid = (int)( fi_endtof / ( 3.141593 / 24 ) );
684 if ( Tofid > 47 ) Tofid = Tofid - 48;
685 tofid_emc[i] = Tofid;
686 if ( module[i] == 2 ) idmatch_emc[2][Tofid] = 1;
687 if ( module[i] == 0 ) idmatch_emc[0][Tofid] = 1;
688 }
689 if ( useEtofMRPC )
690 {
691 int Tofid = (int)( fi_endtof / ( 3.141593 / 18 ) );
692 if ( Tofid > 35 ) Tofid = Tofid - 36;
693 etfid_emc[i] = Tofid;
694 if ( module[i] == 2 ) idmatch_etf_emc[2][Tofid] = 1;
695 if ( module[i] == 0 ) idmatch_etf_emc[0][Tofid] = 1;
696 }
697 }
698 }
699
700 ntrk = 0;
701 if ( ntot > 0 )
702 {
703 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
704 for ( int idfztrk = 0; iter_trk != newtrkCol->end(); iter_trk++, idfztrk++ )
705 {
706 double mdcftrk[5];
707 mdcftrk[0] = ( *iter_trk )->helix( 0 );
708 mdcftrk[1] = ( *iter_trk )->helix( 1 );
709 mdcftrk[2] = -( *iter_trk )->helix( 2 );
710 mdcftrk[3] = ( *iter_trk )->helix( 3 );
711 mdcftrk[4] = ( *iter_trk )->helix( 4 );
712
713 if ( optCosmic == 0 )
714 {
715 Emc_helix EmcHit;
716 // EMC acceptance
717 EmcHit.pathlCut( Estparam.pathlCut() );
718 // Set max. path length(cm)
719 // MDC --> EMC match
720
721 if ( EmcHit.Emc_Get( sqrt( RCEMC2 ), idfztrk, mdcftrk ) > 0 )
722 {
723 double z_emc = EmcHit.Z_emc;
724 double thetaemc_ext = EmcHit.theta_emc;
725 double phiemc_ext = EmcHit.phi_emc;
726
727 double kappa = ( *iter_trk )->helix( 2 );
728 double tanl = ( *iter_trk )->helix( 4 );
729 double _momentum = 1. / fabs( kappa ) * sqrt( 1. + tanl * tanl );
730 for ( int t = 0; t < emctrk; t++ )
731 {
732 if ( ( thetaemc_ext >= ( thetaemc_rec[t] - 0.1 ) ) &&
733 ( thetaemc_ext <= ( thetaemc_rec[t] + 0.1 ) ) &&
734 ( phiemc_ext >= ( phiemc_rec[t] - 0.1 ) ) &&
735 ( phiemc_ext <= ( phiemc_rec[t] + 0.1 ) ) )
736 {
737 if ( ( energy_rec[t] ) >= ( 0.85 * _momentum ) ) particleId[idfztrk] = 1;
738 }
739 }
740 }
741 // get dE/dx results
742 if ( particleId[idfztrk] != 1 )
743 {
744
745 SmartDataPtr<RecMdcDedxCol> newdedxCol( eventSvc(), "/Event/Recon/RecMdcDedxCol" );
746 if ( !newdedxCol || newdedxCol->size() == 0 )
747 { log << MSG::WARNING << "Could not find RecMdcDedxCol" << endmsg; }
748 else
749 {
750 RecMdcDedxCol::iterator it_dedx = newdedxCol->begin();
751 for ( int npid = 0; it_dedx != newdedxCol->end(); it_dedx++, npid++ )
752 {
753 log << MSG::INFO << "retrieved MDC dE/dx: "
754 << "dEdx Id: " << ( *it_dedx )->trackId()
755 << " particle Id: " << ( *it_dedx )->particleType() << endmsg;
756 if ( ( *it_dedx )->particleType() == proton ) { particleId[npid] = 5; }
757 if ( m_debug == 4 ) cout << "dedx pid: " << particleId[npid] << endl;
758 }
759 }
760 }
761 } // if(optCosmic==0)
762
763 idtof = -100;
764 idetf = -100;
765 TofFz_helix TofHit;
766
767 // TOF acceptance
768 TofHit.pathlCut( Estparam.pathlCut() );
769 // Set max. path length(cm)
770
771 TofHit.ztofCut( Estparam.ztofCutmin(),
772 Estparam.ztofCutmax() ); // Set Ztof cut window (cm)
773 // MDC --> TOF match
774 int tofpart = TofHit.TofFz_Get( sqrt( RCTOF2 ), idfztrk, mdcftrk );
775 if ( tofpart < 0 ) continue;
776 // if (TofHit.TofFz_Get(sqrt(RCTOF2),idfztrk,mdcftrk) < 0) continue;
777
778 bool useBarrelScin = ( tofpart == 1 ); // barrel
779 bool useEndcapScin = ( ( tofpart == 0 || tofpart == 2 ) &&
780 useEtofScin ); // endcap for 96Scin and 92Scin+2MRPC;
781 bool useEndcapMRPC = ( ( tofpart == 0 || tofpart == 2 ) &&
782 useEtofMRPC ); // endcap for 72MRPC and 92Scin+2MRPC
783
784 if ( useBarrelScin || useEndcapScin )
785 {
786 idtof = TofHit.Tofid;
787 tofid_helix[idfztrk] = TofHit.Tofid;
788 }
789 if ( useEndcapMRPC )
790 {
791 idetf = TofHit.Etfid;
792 etfid_helix[idfztrk] = TofHit.Etfid;
793 }
794
795 log << MSG::INFO << "helix to tof hit part: " << tofpart << " tof id: " << idtof
796 << " etf id:" << idetf << endmsg;
797 if ( m_debug == 4 )
798 cout << "helix to tof hit part, Id: " << tofpart << " , " << idtof << endl;
799 if ( ( useBarrelScin && idtof >= 0 && idtof <= 87 ) ||
800 ( useEndcapScin && idtof >= 0 && idtof <= 47 ) ||
801 ( useEndcapMRPC && idetf >= 0 && idetf <= 35 ) )
802 { // barrel part: idtof:0~87; endcap part: idtof:0~47 (scintillator),idetf:0~35 (mrpc)
803
804 double abmom = 0.0;
805 if ( useEndcapMRPC )
806 {
807 idetfmatch[tofpart][idetf] = 1;
808 helpathetf[idetf] = TofHit.Path_etf;
809 helz[idetf] = TofHit.Z_etf;
810 abmom = 1. / fabs( TofHit.Kappa ) * sqrt( 1. + TofHit.Tanl * TofHit.Tanl );
811 if ( abmom < 0.1 ) continue;
812 abmom2etf[idetf] = abmom * abmom;
813 r_endetf[idfztrk] = TofHit.r_etf;
814 helzetf[idfztrk] = helz[idetf];
815 }
816
817 if ( useBarrelScin || useEndcapScin )
818 {
819 idmatch[tofpart][idtof] = 1;
820 helpath[idtof] = TofHit.Pathl;
821 helz[idtof] = TofHit.Z_tof;
822 abmom = 1. / fabs( TofHit.Kappa ) * sqrt( 1. + TofHit.Tanl * TofHit.Tanl );
823 if ( abmom < 0.1 ) continue;
824 abmom2[idtof] = abmom * abmom;
825 r_endtof[idfztrk] = TofHit.r_endtof;
826 helztof[idfztrk] = helz[idtof];
827 }
828
829 if ( m_debug == 4 )
830 {
831 cout << "Scintillator info trk number=" << idfztrk << " tofpart=" << tofpart
832 << " idtof=" << idtof << " helpath=" << helpath[idtof]
833 << " helz=" << helz[idtof] << " abmom=" << abmom2[idtof]
834 << " r=" << r_endtof[idfztrk] << " helztof=" << helz[idtof] << endl;
835 cout << "MRPC info trk number=" << idfztrk << " tofpart=" << tofpart
836 << " idetf=" << idetf << " helpath=" << helpathetf[idetf]
837 << " helz=" << helzetf[idetf] << " abmom=" << abmom2etf[idetf]
838 << " r=" << r_endetf[idfztrk] << " helztof=" << helzetf[idetf] << endl;
839 }
840
841 double vel = 1.0e-6;
842 if ( optCosmic == 0 )
843 {
844
845 if ( useEndcapMRPC )
846 {
847 if ( particleId[idfztrk] == 1 )
848 {
849 tetf[idfztrk] =
850 sqrt( ELMAS2 / abmom2etf[idetf] + 1 ) * helpathetf[idetf] / VLIGHT;
851 vel = VLIGHT / sqrt( ELMAS2 / abmom2etf[idetf] + 1 );
852 }
853 else if ( particleId[idfztrk] == 5 )
854 {
855 tetf[idfztrk] =
856 sqrt( PROTONMAS2 / abmom2etf[idetf] + 1 ) * helpathetf[idetf] / VLIGHT;
857 vel = VLIGHT / sqrt( PROTONMAS2 / abmom2etf[idtof] + 1 );
858 }
859 else
860 {
861 tetf[idfztrk] =
862 sqrt( PIMAS2 / abmom2etf[idetf] + 1 ) * helpathetf[idetf] / VLIGHT;
863 vel = VLIGHT / sqrt( PIMAS2 / abmom2etf[idtof] + 1 );
864 }
865 }
866
867 if ( useBarrelScin || useEndcapScin )
868 {
869 if ( particleId[idfztrk] == 1 )
870 {
871 ttof[idfztrk] = sqrt( ELMAS2 / abmom2[idtof] + 1 ) * helpath[idtof] / VLIGHT;
872 vel = VLIGHT / sqrt( ELMAS2 / abmom2[idtof] + 1 );
873 }
874 else if ( particleId[idfztrk] == 5 )
875 {
876 ttof[idfztrk] = sqrt( PROTONMAS2 / abmom2[idtof] + 1 ) * helpath[idtof] / VLIGHT;
877 vel = VLIGHT / sqrt( PROTONMAS2 / abmom2[idtof] + 1 );
878 }
879 else
880 {
881 ttof[idfztrk] = sqrt( PIMAS2 / abmom2[idtof] + 1 ) * helpath[idtof] / VLIGHT;
882 vel = VLIGHT / sqrt( PIMAS2 / abmom2[idtof] + 1 );
883 }
884 }
885 }
886 else
887 {
888
889 if ( useEndcapMRPC )
890 {
891 tetf[idfztrk] = sqrt( MUMAS2 / abmom2etf[idetf] + 1 ) * helpathetf[idetf] / VLIGHT;
892 vel = VLIGHT / sqrt( MUMAS2 / abmom2etf[idtof] + 1 );
893 }
894
895 if ( useBarrelScin || useEndcapMRPC )
896 {
897 ttof[idfztrk] = sqrt( MUMAS2 / abmom2[idtof] + 1 ) * helpath[idtof] / VLIGHT;
898 vel = VLIGHT / sqrt( MUMAS2 / abmom2[idtof] + 1 );
899 }
900 }
901
902 if ( m_ntupleflag && m_tuple2 )
903 {
904 g_vel[idfztrk] = vel;
905 g_abmom[idfztrk] = abmom;
906 if ( useEndcapMRPC ) { g_ttof[idfztrk] = tetf[idfztrk]; }
907 if ( useBarrelScin || useEndcapScin ) { g_ttof[idfztrk] = ttof[idfztrk]; }
908 g_pid[idfztrk] = particleId[idfztrk];
909 }
910 }
911 ntrk++;
912 }
913 if ( m_ntupleflag && m_tuple2 ) g_ntrk = ntrk;
914 }
915
916 // C a l c u l a t e E v t i m e
917 // Retrieve TofMCHit truth
918 if ( runNo < 0 )
919 {
920 SmartDataPtr<TofMcHitCol> tofmcHitCol( eventSvc(), "/Event/MC/TofMcHitCol" );
921 if ( !tofmcHitCol )
922 {
923 log << MSG::ERROR << "Could not find McParticle" << endmsg;
924 // return StatusCode::FAILURE;
925 }
926 else
927 {
928 TofMcHitCol::iterator iter_mctof = tofmcHitCol->begin();
929
930 for ( ; iter_mctof != tofmcHitCol->end(); iter_mctof++, digiId++ )
931 {
932 log << MSG::INFO << " TofMcHit Flight Time = " << ( *iter_mctof )->getFlightTime()
933 << " zPosition = " << ( ( *iter_mctof )->getPositionZ() ) / 10
934 << " xPosition = " << ( ( *iter_mctof )->getPositionX() ) / 10
935 << " yPosition = " << ( ( *iter_mctof )->getPositionY() ) / 10 << endmsg;
936 }
937 }
938 }
939
940 ////choose cosmic
941 TofDataVector tofDigiVec = m_rawDataProviderSvc->tofDataVectorEstime();
942 // if(tofDigiVec.size()==0) return StatusCode::SUCCESS;
943 for ( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end();
944 iter2++ )
945 {
946 int barrelid;
947 int layerid;
948 int tofid;
949 int strip;
950
951 if ( !( ( *iter2 )->is_mrpc() ) )
952 { // Scintillator
953 if ( ( *iter2 )->barrel() )
954 { // Scintillator Barrel
955 barrelid = 1;
956 tofid = ( *iter2 )->tofId();
957 layerid = ( *iter2 )->layer();
958 if ( layerid == 1 ) tofid = tofid - 88;
959 if ( ( ( *iter2 )->quality() & 0x5 ) == 0x5 && ( *iter2 )->times() == 1 )
960 {
961 double ftdc = ( *iter2 )->tdc1();
962 double btdc = ( *iter2 )->tdc2();
963 mean_tdc_btof[layerid][tofid] = ( ftdc + btdc ) / 2;
964 }
965 else if ( ( ( *iter2 )->quality() & 0x5 ) == 0x5 && ( *iter2 )->times() > 1 )
966 {
967 double ftdc = ( *iter2 )->tdc1();
968 double btdc = ( *iter2 )->tdc2();
969 mean_tdc_btof[layerid][tofid] = ( ftdc + btdc ) / 2;
970 }
971 }
972 else
973 { // Scintillator Endcap
974 tofid = ( *iter2 )->tofId();
975 if ( tofid < 48 ) barrelid = 0;
976 if ( tofid > 47 ) barrelid = 2;
977 if ( barrelid == 2 ) tofid = tofid - 48;
978
979 if ( ( *iter2 )->times() == 1 )
980 {
981 double ftdc = ( *iter2 )->tdc();
982 mean_tdc_etof[barrelid][tofid] = ftdc;
983 }
984 else if ( ( ( *iter2 )->times() > 1 ) && ( ( *iter2 )->times() < 5 ) )
985 {
986 double ftdc = ( *iter2 )->tdc();
987 mean_tdc_etof[barrelid][tofid] = ftdc;
988 }
989 }
990 }
991 else
992 { // MRPC Endcap
993 tofid = ( *iter2 )->tofId();
994 strip = ( *iter2 )->strip();
995 if ( tofid < 36 ) barrelid = 0;
996 if ( tofid > 35 ) barrelid = 2;
997 if ( barrelid == 2 ) tofid = tofid - 36;
998 if ( ( ( *iter2 )->quality() & 0x5 ) == 0x5 && ( *iter2 )->times() == 1 )
999 {
1000 double ftdc = ( *iter2 )->tdc1();
1001 double btdc = ( *iter2 )->tdc2();
1002 mean_tdc_etf[barrelid][tofid][strip] = ( ftdc + btdc ) / 2;
1003 }
1004 else if ( ( ( *iter2 )->quality() & 0x5 ) == 0x5 && ( *iter2 )->times() > 1 )
1005 {
1006 double ftdc = ( *iter2 )->tdc1();
1007 double btdc = ( *iter2 )->tdc2();
1008 mean_tdc_etf[barrelid][tofid][strip] = ( ftdc + btdc ) / 2;
1009 }
1010 }
1011 }
1012
1013 double difftof_b = 100, difftof_e = 100;
1014 int tofid1 = tofid_emc[0];
1015 int tofid2 = tofid_emc[1];
1016 if ( module[0] == 1 && module[1] == 1 )
1017 {
1018 for ( int i = 0; i < 2; i++ )
1019 {
1020 for ( int m = 0; m < 2; m++ )
1021 {
1022 for ( int j = -2; j < 3; j++ )
1023 {
1024 for ( int k = -2; k < 3; k++ )
1025 {
1026 int p = tofid1 + j;
1027 int q = tofid2 + k;
1028 if ( p < 0 ) p = p + 88;
1029 if ( p > 87 ) p = p - 88;
1030 if ( q < 0 ) q = q + 88;
1031 if ( q > 87 ) q = q + 88;
1032 if ( mean_tdc_btof[i][p] == 0 || mean_tdc_btof[m][q] == 0 ) continue;
1033 double difftof_b_temp = mean_tdc_btof[i][p] - mean_tdc_btof[m][q];
1034 if ( abs( difftof_b_temp ) < abs( difftof_b ) ) difftof_b = difftof_b_temp;
1035 if ( m_ntupleflag && m_tuple2 ) g_difftof_b = difftof_b;
1036 }
1037 }
1038 }
1039 }
1040 }
1041
1042 if ( useEtofMRPC )
1043 {
1044 if ( module[0] != 1 && module[1] != 1 )
1045 {
1046 tofid1 = etfid_emc[0];
1047 tofid2 = etfid_emc[1];
1048 for ( int i = -1; i < 2; i++ )
1049 {
1050 for ( int j = -1; j < 2; j++ )
1051 {
1052 int m = tofid1 + i;
1053 int n = tofid2 + j;
1054 if ( m < 0 ) m = 36 + m;
1055 if ( m > 35 ) m = m - 36;
1056 if ( n < 0 ) n = 36 + n;
1057 if ( n > 35 ) n = n - 36;
1058 if ( mean_tdc_etf[0][m] && mean_tdc_etf[2][n] )
1059 {
1060 double difftof_e_temp = mean_tdc_etf[0][m] - mean_tdc_etf[2][n];
1061 if ( abs( difftof_e_temp ) < abs( difftof_e ) ) difftof_e = difftof_e_temp;
1062 if ( m_ntupleflag && m_tuple2 ) g_difftof_e = difftof_e;
1063 }
1064 }
1065 }
1066 }
1067 }
1068
1069 if ( useEtofScin )
1070 {
1071 if ( module[0] != 1 && module[1] != 1 )
1072 {
1073 for ( int i = -1; i < 2; i++ )
1074 {
1075 for ( int j = -1; j < 2; j++ )
1076 {
1077 int m = tofid1 + i;
1078 int n = tofid2 + j;
1079 if ( m < 0 ) m = 48 + m;
1080 if ( m > 47 ) m = m - 48;
1081 if ( n < 0 ) n = 48 + n;
1082 if ( n > 47 ) n = n - 48;
1083 if ( mean_tdc_etof[0][m] && mean_tdc_etof[2][n] )
1084 {
1085 double difftof_e_temp = mean_tdc_etof[0][m] - mean_tdc_etof[2][n];
1086 if ( abs( difftof_e_temp ) < abs( difftof_e ) ) difftof_e = difftof_e_temp;
1087 if ( m_ntupleflag && m_tuple2 ) g_difftof_e = difftof_e;
1088 }
1089 }
1090 }
1091 }
1092 }
1093
1094 if ( m_optCosmic == 1 ) optCosmic = 1;
1095 else optCosmic = 0;
1096
1097 digiId = 0;
1098 unsigned int tofid;
1099 unsigned int barrelid;
1100 unsigned int layerid;
1101 t0forward_add = 0;
1102 t0backward_add = 0;
1103 TofDataVector::iterator iter2 = tofDigiVec.begin();
1104 for ( ; iter2 != tofDigiVec.end(); iter2++, digiId++ )
1105 {
1106 log << MSG::INFO << "TOF digit No: " << digiId << endmsg;
1107 barrelid = ( *iter2 )->barrel();
1108 if ( ( *iter2 )->barrel() == 0 ) continue;
1109 if ( ( ( *iter2 )->quality() & 0x5 ) == 0x5 && ( *iter2 )->times() == 1 )
1110 { // tof_flag=1
1111 tofid = ( *iter2 )->tofId();
1112 layerid = ( *iter2 )->layer();
1113 if ( layerid == 1 ) tofid = tofid - 88;
1114 log << MSG::INFO << " TofId = " << tofid << " barrelid = " << barrelid
1115 << " layerid = " << layerid << " ForwordADC = " << ( *iter2 )->adc1()
1116 << " ForwordTDC = " << ( *iter2 )->tdc1()
1117 << " BackwordADC = " << ( *iter2 )->adc2()
1118 << " BackwordTDC = " << ( *iter2 )->tdc2() << endmsg;
1119 // forb=0/1 for forward(z>0, east) and backward(z<0, west)
1120 double ftdc = ( *iter2 )->tdc1();
1121 double btdc = ( *iter2 )->tdc2();
1122 if ( m_debug == 4 )
1123 cout << "barrel 1 ::layer, barrel, tofid, ftdc, btdc: " << layerid << " , " << barrelid
1124 << " , " << tofid << " , " << ftdc << " , " << btdc << endl;
1125 double fadc = ( *iter2 )->adc1();
1126 double badc = ( *iter2 )->adc2();
1127 int idptof = ( ( tofid - 1 ) == -1 ) ? 87 : tofid - 1;
1128 int idntof = ( ( tofid + 1 ) == 88 ) ? 0 : tofid + 1;
1129 double ztof = fabs( ( ftdc - btdc ) / 2 ) * 17.7, ztof2 = ztof * ztof;
1130 double mean_tdc = 0.5 * ( btdc + ftdc );
1131 double cor_path = sqrt( RCTOF2 + ztof2 ) / VLIGHT;
1132
1133 if ( idmatch[barrelid][tofid] == 1 || idmatch[barrelid][idptof] == 1 ||
1134 idmatch[barrelid][idntof] == 1 )
1135 {
1136 for ( int i = 0; i < ntot; i++ )
1137 {
1138 if ( ttof[i] != 0 && ftdc > 0 )
1139 {
1140 if ( ( tofid_helix[i] == tofid ) || ( tofid_helix[i] == idntof ) ||
1141 ( tofid_helix[i] == idptof ) )
1142 {
1143 if ( barrelid == 1 && helztof[i] < 117 && helztof[i] > -117 )
1144 {
1145 if ( optCosmic && tofid < 44 )
1146 {
1147 backevtime = -ttof[i] + ( 115 + helztof[i] ) * 0.0566 + 12.;
1148 forevtime = -ttof[i] + ( 115 - helztof[i] ) * 0.0566 + 12.;
1149 meantevup[ntofup] = ( backevtime + forevtime ) / 2;
1150 ntofup++;
1151 }
1152 else
1153 {
1154 backevtime = ttof[i] + ( 115 + helztof[i] ) * 0.0566 + 12.;
1155 forevtime = ttof[i] + ( 115 - helztof[i] ) * 0.0566 + 12.;
1156 meantevdown[ntofdown] = ( backevtime + forevtime ) / 2;
1157 ntofdown++;
1158 }
1159 if ( ( *iter2 )->adc1() < 0.0 || ( *iter2 )->adc2() < 0.0 || m_userawtime )
1160 {
1161 t0forward_trk = ftdc - forevtime;
1162 t0backward_trk = btdc - backevtime;
1163 }
1164 else
1165 {
1166 t0forward_trk = tofCaliSvc->BTime1( ( *iter2 )->adc1(), ( *iter2 )->tdc1(),
1167 helztof[i], ( *iter2 )->tofId() ) -
1168 ttof[i];
1169 t0backward_trk = tofCaliSvc->BTime2( ( *iter2 )->adc2(), ( *iter2 )->tdc2(),
1170 helztof[i], ( *iter2 )->tofId() ) -
1171 ttof[i];
1172 if ( optCosmic && tofid < 44 )
1173 {
1174 t0forward_trk = tofCaliSvc->BTime1( ( *iter2 )->adc1(), ( *iter2 )->tdc1(),
1175 helztof[i], ( *iter2 )->tofId() ) +
1176 ttof[i];
1177 t0backward_trk =
1178 tofCaliSvc->BTime2( ( *iter2 )->adc2(), ( *iter2 )->tdc2(), helztof[i],
1179 ( *iter2 )->tofId() ) +
1180 ttof[i];
1181 }
1182 }
1183
1184 if ( t0forward_trk < -3 || t0backward_trk < -3 ||
1185 fabs( t0forward_trk - t0backward_trk ) > 15.0 )
1186 continue;
1187 if ( !m_TofOpt && nmatch_barrel != 0 &&
1188 fabs( ( t0forward_trk + t0backward_trk ) / 2 -
1189 ( t0backward_add + t0forward_add ) / 2 / nmatch_barrel ) > 11 )
1190 continue;
1191 if ( m_debug == 4 )
1192 cout << " t0forward_trk, t0backward_trk: " << t0forward_trk << " , "
1193 << t0backward_trk << endl;
1194 if ( m_ntupleflag && m_tuple2 )
1195 {
1196 g_t0for[nmatch1] = t0forward_trk;
1197 g_t0back[nmatch2] = t0backward_trk;
1198 g_meantdc = ( ftdc + btdc ) / 2;
1199 if ( tofid < 44 ) g_ntofup1++;
1200 else g_ntofdown1++;
1201 }
1202 t0forward_add += t0forward_trk;
1203 t0backward_add += t0backward_trk;
1204 if ( nmatch > 499 ) break;
1205 meantev[nmatch] = ( backevtime + forevtime ) / 2;
1206 Tof_t0Arr[nmatch] = ( t0forward_trk + t0backward_trk ) / 2.0;
1207 nmatch++;
1208 nmatch1 = nmatch1 + 1;
1209 nmatch2 = nmatch2 + 1;
1210 nmatch_barrel++;
1211 }
1212 }
1213 }
1214 }
1215 }
1216 }
1217 }
1218
1219 if ( nmatch_barrel != 0 )
1220 {
1221 if ( m_ntupleflag && m_tuple2 )
1222 { g_t0barrelTof = ( t0forward_add / nmatch_barrel + t0backward_add / nmatch_barrel ) / 2; }
1223 tof_flag = 1;
1224 t_quality = 1;
1225 }
1226
1227 if ( nmatch_barrel == 0 )
1228 {
1229 digiId = 0;
1230 t0forward_add = 0;
1231 t0backward_add = 0;
1232 for ( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end();
1233 iter2++, digiId++ )
1234 {
1235 log << MSG::INFO << "TOF digit No: " << digiId << endmsg;
1236 barrelid = ( *iter2 )->barrel();
1237 if ( ( *iter2 )->barrel() == 0 ) continue;
1238 if ( ( ( *iter2 )->quality() & 0x5 ) == 0x5 && ( *iter2 )->times() > 1 )
1239 { // tof_flag=2
1240 tofid = ( *iter2 )->tofId();
1241 layerid = ( *iter2 )->layer();
1242 if ( layerid == 1 ) tofid = tofid - 88;
1243 log << MSG::INFO << " TofId = " << tofid << " barrelid = " << barrelid
1244 << " layerid = " << layerid << " ForwordADC = " << ( *iter2 )->adc1()
1245 << " ForwordTDC = " << ( *iter2 )->tdc1() << endmsg;
1246 double ftdc = ( *iter2 )->tdc1();
1247 double btdc = ( *iter2 )->tdc2();
1248 double fadc = ( *iter2 )->adc1();
1249 double badc = ( *iter2 )->adc2();
1250 if ( m_debug == 4 )
1251 cout << "barrel 2 ::layer, barrel, tofid, ftdc, btdc: " << layerid << " , "
1252 << barrelid << " , " << tofid << " , " << ftdc << " , " << btdc << endl;
1253 int idptof = ( ( tofid - 1 ) == -1 ) ? 87 : tofid - 1;
1254 int idntof = ( ( tofid + 1 ) == 88 ) ? 0 : tofid + 1;
1255 if ( idmatch[barrelid][tofid] == 1 || idmatch[barrelid][idptof] == 1 ||
1256 idmatch[barrelid][idntof] == 1 )
1257 {
1258 for ( int i = 0; i < ntot; i++ )
1259 {
1260 if ( ttof[i] != 0 && ftdc > 0 )
1261 {
1262 if ( tofid_helix[i] == tofid || ( tofid_helix[i] == idptof ) ||
1263 ( tofid_helix[i] == idntof ) )
1264 {
1265 if ( barrelid == 1 && helztof[i] < 117 && helztof[i] > -117 )
1266 {
1267 if ( optCosmic && tofid < 44 )
1268 {
1269 backevtime = -ttof[i] + ( 115 + helztof[i] ) * 0.0566 + 12.;
1270 forevtime = -ttof[i] + ( 115 - helztof[i] ) * 0.0566 + 12.;
1271 meantevup[ntofup] = ( backevtime + forevtime ) / 2;
1272 ntofup++;
1273 }
1274 else
1275 {
1276 backevtime = ttof[i] + ( 115 + helztof[i] ) * 0.0566 + 12.;
1277 forevtime = ttof[i] + ( 115 - helztof[i] ) * 0.0566 + 12.;
1278 meantevdown[ntofdown] = ( backevtime + forevtime ) / 2;
1279 ntofdown++;
1280 }
1281 if ( ( *iter2 )->adc1() < 0.0 || ( *iter2 )->adc2() < 0.0 || m_userawtime )
1282 {
1283 t0forward_trk = ftdc - forevtime;
1284 t0backward_trk = btdc - backevtime;
1285 }
1286 else
1287 {
1288 t0forward_trk = tofCaliSvc->BTime1( ( *iter2 )->adc1(), ( *iter2 )->tdc1(),
1289 helztof[i], ( *iter2 )->tofId() ) -
1290 ttof[i];
1291 t0backward_trk =
1292 tofCaliSvc->BTime2( ( *iter2 )->adc2(), ( *iter2 )->tdc2(), helztof[i],
1293 ( *iter2 )->tofId() ) -
1294 ttof[i];
1295 if ( optCosmic && tofid < 44 )
1296 {
1297 t0forward_trk =
1298 tofCaliSvc->BTime1( ( *iter2 )->adc1(), ( *iter2 )->tdc1(),
1299 helztof[i], ( *iter2 )->tofId() ) +
1300 ttof[i];
1301 t0backward_trk =
1302 tofCaliSvc->BTime2( ( *iter2 )->adc2(), ( *iter2 )->tdc2(),
1303 helztof[i], ( *iter2 )->tofId() ) +
1304 ttof[i];
1305 }
1306 }
1307
1308 if ( t0forward_trk < -3 || t0backward_trk < -3 ||
1309 fabs( t0forward_trk - t0backward_trk ) > 15.0 )
1310 continue;
1311 if ( !m_TofOpt && nmatch_barrel != 0 &&
1312 fabs( ( t0forward_trk + t0backward_trk ) / 2 -
1313 ( t0backward_add + t0forward_add ) / 2 / nmatch_barrel ) > 11 )
1314 continue;
1315 if ( m_debug == 4 )
1316 cout << "t0forward_trk, t0backward_trk: " << t0forward_trk << " , "
1317 << t0backward_trk << endl;
1318 if ( m_ntupleflag && m_tuple2 )
1319 {
1320 g_t0for[nmatch1] = t0forward_trk;
1321 g_t0back[nmatch2] = t0backward_trk;
1322 g_meantdc = ( ftdc + btdc ) / 2;
1323 if ( tofid < 44 ) g_ntofup1++;
1324 else g_ntofdown1++;
1325 }
1326 t0forward_add += t0forward_trk;
1327 t0backward_add += t0backward_trk;
1328 if ( nmatch > 499 ) break;
1329 meantev[nmatch] = ( backevtime + forevtime ) / 2;
1330 Tof_t0Arr[nmatch] = ( t0forward_trk + t0backward_trk ) / 2.0;
1331 nmatch++;
1332 nmatch1 = nmatch1 + 1;
1333 nmatch2 = nmatch2 + 1;
1334 nmatch_barrel++;
1335 }
1336 }
1337 }
1338 }
1339 }
1340 }
1341 } // close 2nd for loop -> only barrel part
1342 if ( nmatch_barrel ) tof_flag = 2;
1343 }
1344
1345 if ( ntot == 0 || nmatch_barrel == 0 )
1346 {
1347 for ( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end();
1348 iter2++, digiId++ )
1349 {
1350 log << MSG::INFO << "TOF digit No: " << digiId << endmsg;
1351 barrelid = ( *iter2 )->barrel();
1352 if ( ( *iter2 )->barrel() == 0 ) continue;
1353 if ( ( ( *iter2 )->quality() & 0x5 ) == 0x5 && ( *iter2 )->times() == 1 )
1354 { // tof_flag=3
1355 tofid = ( *iter2 )->tofId();
1356 layerid = ( *iter2 )->layer();
1357 if ( layerid == 1 ) tofid = tofid - 88;
1358 log << MSG::INFO << " TofId = " << tofid << " barrelid = " << barrelid
1359 << " layerid = " << layerid << " ForwordADC = " << ( *iter2 )->adc1()
1360 << " ForwordTDC = " << ( *iter2 )->tdc1() << endmsg;
1361 double ftdc = ( *iter2 )->tdc1();
1362 double btdc = ( *iter2 )->tdc2();
1363 double fadc = ( *iter2 )->adc1();
1364 double badc = ( *iter2 )->adc2();
1365 if ( m_debug == 4 )
1366 cout << "barrel 3 ::layer, barrel, tofid, ftdc, btdc: " << layerid << " , "
1367 << barrelid << " , " << tofid << " , " << ftdc << " , " << btdc << endl;
1368 int idptof = ( ( tofid - 1 ) == -1 ) ? 87 : tofid - 1;
1369 int idntof = ( ( tofid + 1 ) == 88 ) ? 0 : tofid + 1;
1370 for ( int i = 0; i < 2; i++ )
1371 {
1372 if ( tofid_emc[i] == tofid || tofid_emc[i] == idptof || tofid_emc[i] == idntof )
1373 {
1374 if ( zemc_rec[0] || zemc_rec[1] )
1375 {
1376 if ( tofid == tofid_emc[i] || tofid_emc[i] == idntof || tofid_emc[i] == idptof )
1377 {
1378 if ( ftdc > 2000. || module[i] != 1 ) continue;
1379 ttof[i] = sqrt( ELMAS2 / ( m_ebeam * m_ebeam ) + 1 ) *
1380 sqrt( xemc_rec[i] * xemc_rec[i] + yemc_rec[i] * yemc_rec[i] +
1381 zemc_rec[i] * zemc_rec[i] ) /
1382 VLIGHT;
1383 if ( optCosmic == 1 && tofid < 44 )
1384 {
1385 backevtime = -ttof[i] + ( 115 + zemc_rec[i] ) * 0.0566 + 12.;
1386 forevtime = -ttof[i] + ( 115 - zemc_rec[i] ) * 0.0566 + 12.;
1387 meantevup[ntofup] = ( backevtime + forevtime ) / 2;
1388 ntofup++;
1389 }
1390 else
1391 {
1392 backevtime = ttof[i] + ( 115 + zemc_rec[i] ) * 0.0566 + 12.;
1393 forevtime = ttof[i] + ( 115 - zemc_rec[i] ) * 0.0566 + 12.;
1394 meantevdown[ntofdown] = ( backevtime + forevtime ) / 2;
1395 ntofdown++;
1396 }
1397 if ( ( *iter2 )->adc1() < 0.0 || ( *iter2 )->adc2() < 0.0 || m_userawtime )
1398 {
1399 t0forward_trk = ftdc - forevtime;
1400 t0backward_trk = btdc - backevtime;
1401 }
1402 else
1403 {
1404 t0forward_trk = tofCaliSvc->BTime1( ( *iter2 )->adc1(), ( *iter2 )->tdc1(),
1405 helztof[i], ( *iter2 )->tofId() ) -
1406 ttof[i];
1407 t0backward_trk = tofCaliSvc->BTime2( ( *iter2 )->adc2(), ( *iter2 )->tdc2(),
1408 helztof[i], ( *iter2 )->tofId() ) -
1409 ttof[i];
1410 if ( optCosmic && tofid < 44 )
1411 {
1412 t0forward_trk = tofCaliSvc->BTime1( ( *iter2 )->adc1(), ( *iter2 )->tdc1(),
1413 helztof[i], ( *iter2 )->tofId() ) +
1414 ttof[i];
1415 t0backward_trk =
1416 tofCaliSvc->BTime2( ( *iter2 )->adc2(), ( *iter2 )->tdc2(), helztof[i],
1417 ( *iter2 )->tofId() ) +
1418 ttof[i];
1419 }
1420 }
1421
1422 if ( t0forward_trk < -1 || t0backward_trk < -1 ||
1423 fabs( t0forward_trk - t0backward_trk ) > 15.0 )
1424 continue;
1425 if ( !m_TofOpt && nmatch_barrel != 0 &&
1426 fabs( ( t0forward_trk + t0backward_trk ) / 2 -
1427 ( t0backward_add + t0forward_add ) / 2 / nmatch_barrel ) > 11 )
1428 continue;
1429 if ( m_debug == 4 )
1430 cout << "t0forward_trk, t0backward_trk: " << t0forward_trk << " , "
1431 << t0backward_trk << endl;
1432 t0forward_add += t0forward_trk;
1433 t0backward_add += t0backward_trk;
1434 if ( nmatch > 499 ) break;
1435 meantev[nmatch] = ( backevtime + forevtime ) / 2;
1436 Tof_t0Arr[nmatch] = ( t0forward_trk + t0backward_trk ) / 2.0;
1437 nmatch++;
1438 nmatch_barrel++;
1439 emcflag1 = 1;
1440 }
1441 }
1442 }
1443 }
1444 }
1445 } // 3rd for loop only barrel part
1446 if ( nmatch_barrel ) tof_flag = 3;
1447 }
1448
1449 if ( nmatch_barrel != 0 )
1450 { // only matched barrel tracks
1451 t0forward = t0forward_add / nmatch_barrel;
1452 t0backward = t0backward_add / nmatch_barrel;
1453 if ( optCosmic == 0 )
1454 {
1455 if ( m_TofOpt )
1456 {
1457 t_Est = EST_Trimmer( Opt_new( Tof_t0Arr, nmatch, m_TofOpt_Cut ), tOffset_b,
1458 toffset_rawtime, offset_dt, bunchtime );
1459 }
1460 else
1461 t_Est = EST_Trimmer( ( t0forward + t0backward ) / 2, tOffset_b, toffset_rawtime,
1462 offset_dt, bunchtime );
1463 if ( t_Est < 0 ) t_Est = 0;
1464 if ( tof_flag == 1 ) tEstFlag = 111;
1465 else if ( tof_flag == 2 ) tEstFlag = 121;
1466 else if ( tof_flag == 3 ) tEstFlag = 131;
1467 }
1468 if ( optCosmic )
1469 {
1470 t_Est = ( t0forward + t0backward ) / 2; // 3. yzhang add tOffset_b for barrel cosmic
1471 if ( tof_flag == 1 ) tEstFlag = 211;
1472 else if ( tof_flag == 2 ) tEstFlag = 221;
1473 else if ( tof_flag == 3 ) tEstFlag = 231;
1474 }
1475 if ( m_ntupleflag && m_tuple2 ) g_meant0 = ( t0forward + t0backward ) / 2;
1476 } // close if(nmatch_barrel !=0)
1477
1478 digiId = 0;
1479 t0forward_add = 0;
1480 t0backward_add = 0;
1481
1482 if ( nmatch_barrel == 0 )
1483 {
1484 for ( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end();
1485 iter2++, digiId++ )
1486 {
1487 log << MSG::INFO << "TOF digit No: " << digiId << endmsg;
1488 barrelid = ( *iter2 )->barrel();
1489 if ( ( *iter2 )->barrel() == 0 ) continue;
1490 if ( ( ( *iter2 )->quality() & 0x5 ) == 0x4 )
1491 { // tEstFlag=241
1492 tofid = ( *iter2 )->tofId();
1493 layerid = ( *iter2 )->layer();
1494 if ( layerid == 1 ) tofid = tofid - 88;
1495 log << MSG::INFO << " TofId = " << tofid << " barrelid = " << barrelid
1496 << " layerid = " << layerid << " ForwordADC = " << ( *iter2 )->adc1()
1497 << " ForwordTDC = " << ( *iter2 )->tdc1() << endmsg;
1498 // forb=0/1 for forward(z>0, east) and backward(z<0, west)
1499 double ftdc = ( *iter2 )->tdc1();
1500 double fadc = ( *iter2 )->adc1();
1501 if ( m_debug == 4 )
1502 cout << "barrel 4 ::layer, barrel, tofid, ftdc: " << layerid << " , " << barrelid
1503 << " , " << tofid << " , " << ftdc << endl;
1504 int idptof = ( ( tofid - 1 ) == -1 ) ? 87 : tofid - 1;
1505 int idntof = ( ( tofid + 1 ) == 88 ) ? 0 : tofid + 1;
1506 if ( idmatch[barrelid][tofid] == 1 || idmatch[barrelid][idptof] == 1 ||
1507 idmatch[barrelid][idntof] == 1 )
1508 {
1509 for ( int i = 0; i < ntot; i++ )
1510 {
1511 if ( ttof[i] != 0 && ftdc > 0 )
1512 {
1513 if ( tofid_helix[i] == tofid || ( tofid_helix[i] == idptof ) ||
1514 ( tofid_helix[i] == idntof ) )
1515 {
1516 if ( barrelid == 1 && helztof[i] < 117 && helztof[i] > -117 )
1517 {
1518 if ( optCosmic && tofid < 44 )
1519 {
1520 forevtime = -ttof[i] + ( 115 - helztof[i] ) * 0.0566 + 12.;
1521 meantevup[ntofup] = forevtime;
1522 ntofup++;
1523 }
1524 else
1525 {
1526 forevtime = ttof[i] + ( 115 - helztof[i] ) * 0.0566 + 12.;
1527 meantevdown[ntofdown] = forevtime;
1528 ntofdown++;
1529 }
1530 if ( ( *iter2 )->adc1() < 0.0 || m_userawtime )
1531 { t0forward_trk = ftdc - forevtime; }
1532 else
1533 {
1534 t0forward_trk = tofCaliSvc->BTime1( ( *iter2 )->adc1(), ( *iter2 )->tdc1(),
1535 helztof[i], ( *iter2 )->tofId() ) -
1536 ttof[i];
1537 }
1538
1539 if ( t0forward_trk < -1 ) continue;
1540 if ( !m_TofOpt && nmatch_barrel_1 != 0 &&
1541 fabs( ( t0forward_trk ) - ( t0forward_add ) / nmatch_barrel_1 ) > 11 )
1542 continue;
1543 if ( m_debug == 4 ) cout << "t0forward_trk: " << t0forward_trk << endl;
1544 if ( m_ntupleflag && m_tuple2 )
1545 {
1546 g_t0for[nmatch1] = t0forward_trk;
1547 g_meantdc = ftdc;
1548 if ( tofid < 44 ) g_ntofup1++;
1549 else g_ntofdown1++;
1550 }
1551 t0forward_add += t0forward_trk;
1552 // t0v.push_back(t0forward_trk);
1553 if ( nmatch > 499 ) break;
1554 meantev[nmatch] = forevtime;
1555 Tof_t0Arr[nmatch] = t0forward_trk;
1556 nmatch++;
1557 nmatch1++;
1558 nmatch_barrel_1++;
1559 }
1560 }
1561 }
1562 }
1563 }
1564 }
1565 else if ( ( ( *iter2 )->quality() & 0x5 ) == 0x1 )
1566 { // tEstFlag=241
1567 tofid = ( *iter2 )->tofId();
1568 layerid = ( *iter2 )->layer();
1569 if ( layerid == 1 ) tofid = tofid - 88;
1570 log << MSG::INFO << " TofId = " << tofid << " barrelid = " << barrelid
1571 << " layerid = " << layerid << " BackwordADC = " << ( *iter2 )->adc2()
1572 << " BackwordTDC = " << ( *iter2 )->tdc2() << endmsg;
1573 // forb=0/1 for forward(z>0, east) and backward(z<0, west)
1574 double btdc = ( *iter2 )->tdc2();
1575 if ( m_debug == 4 )
1576 cout << "barrel 5 ::layer, barrel, tofid, btdc: " << layerid << " , " << barrelid
1577 << " , " << tofid << " , " << btdc << endl;
1578 double badc = ( *iter2 )->adc2();
1579 int idptof = ( ( tofid - 1 ) == -1 ) ? 87 : tofid - 1;
1580 int idntof = ( ( tofid + 1 ) == 88 ) ? 0 : tofid + 1;
1581 if ( idmatch[barrelid][tofid] == 1 || idmatch[barrelid][idptof] == 1 ||
1582 idmatch[barrelid][idntof] == 1 )
1583 {
1584 for ( int i = 0; i < ntot; i++ )
1585 {
1586 if ( ttof[i] != 0 && btdc > 0 )
1587 {
1588 if ( ( tofid_helix[i] == tofid ) || ( tofid_helix[i] == idntof ) ||
1589 ( tofid_helix[i] == idptof ) )
1590 {
1591 if ( barrelid == 1 && helztof[i] < 117 && helztof[i] > -117 )
1592 {
1593 if ( optCosmic && tofid < 44 )
1594 {
1595 backevtime = -ttof[i] + ( 115 + helztof[i] ) * 0.0566 + 12.;
1596 meantevup[ntofup] = backevtime;
1597 ntofup++;
1598 }
1599 else
1600 {
1601 backevtime = ttof[i] + ( 115 + helztof[i] ) * 0.0566 + 12.;
1602 meantevdown[ntofdown] = backevtime;
1603 ntofdown++;
1604 }
1605
1606 if ( ( *iter2 )->adc2() < 0.0 || m_userawtime )
1607 { t0backward_trk = btdc - backevtime; }
1608 else
1609 {
1610 t0backward_trk =
1611 tofCaliSvc->BTime2( ( *iter2 )->adc2(), ( *iter2 )->tdc2(), helztof[i],
1612 ( *iter2 )->tofId() ) -
1613 ttof[i];
1614 }
1615
1616 if ( t0backward_trk < -1 ) continue;
1617 if ( !m_TofOpt && nmatch_barrel_2 != 0 &&
1618 fabs( ( t0backward_trk ) - ( t0backward_add ) / nmatch_barrel_2 ) > 11 )
1619 continue;
1620 if ( m_debug == 4 ) cout << "t0backward_trk: " << t0backward_trk << endl;
1621 if ( m_ntupleflag && m_tuple2 )
1622 {
1623 g_t0back[nmatch2] = t0backward_trk;
1624 g_meantdc = btdc;
1625 if ( tofid < 44 ) g_ntofup1++;
1626 else g_ntofdown1++;
1627 }
1628 t0backward_add += t0backward_trk;
1629 if ( nmatch > 499 ) break;
1630 meantev[nmatch] = backevtime;
1631 Tof_t0Arr[nmatch] = t0backward_trk;
1632 nmatch++;
1633 nmatch2++;
1634 nmatch_barrel_2++;
1635 }
1636 }
1637 }
1638 }
1639 }
1640 }
1641 }
1642 }
1643 if ( nmatch_barrel_1 != 0 )
1644 {
1645 t0forward = t0forward_add / nmatch_barrel_1;
1646 if ( optCosmic == 0 )
1647 {
1648 if ( m_TofOpt )
1649 {
1650 t_Est = EST_Trimmer( Opt_new( Tof_t0Arr, nmatch, m_TofOpt_Cut ), tOffset_b,
1651 toffset_rawtime, offset_dt, bunchtime );
1652 }
1653 else t_Est = EST_Trimmer( t0forward, tOffset_b, toffset_rawtime, offset_dt, bunchtime );
1654 if ( t_Est < 0 ) t_Est = 0;
1655 tEstFlag = 141;
1656 }
1657 else
1658 {
1659 t_Est = t0forward; // 5. yzhang add for nmatch_barrel_1 cosmic
1660 tEstFlag = 241;
1661 }
1662 if ( m_ntupleflag && m_tuple2 ) g_meant0 = t0forward;
1663 }
1664 if ( nmatch_barrel_2 != 0 )
1665 {
1666 t0backward = t0backward_add / nmatch_barrel_2;
1667 if ( optCosmic == 0 )
1668 {
1669 if ( m_TofOpt )
1670 {
1671 t_Est = EST_Trimmer( Opt_new( Tof_t0Arr, nmatch, m_TofOpt_Cut ), tOffset_b,
1672 toffset_rawtime, offset_dt, bunchtime );
1673 }
1674 else t_Est = EST_Trimmer( t0backward, tOffset_b, toffset_rawtime, offset_dt, bunchtime );
1675 if ( t_Est < 0 ) t_Est = 0;
1676 tEstFlag = 141;
1677 }
1678 else
1679 {
1680 t_Est = t0backward; // 7. yzhang add tOffset_b for nmatch_barrel_2 cosmic
1681 tEstFlag = 241;
1682 }
1683 if ( m_ntupleflag && m_tuple2 ) g_meant0 = t0backward;
1684 }
1685
1686 digiId = 0;
1687 t0forward_add = 0;
1688 t0backward_add = 0;
1689 if ( nmatch_barrel == 0 )
1690 {
1691 for ( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end();
1692 iter2++, digiId++ )
1693 {
1694 log << MSG::INFO << "TOF digit No: " << digiId << endmsg;
1695 barrelid = ( *iter2 )->barrel();
1696 if ( ( *iter2 )->barrel() != 0 ) continue;
1697 if ( ( *iter2 )->times() != 1 ) continue;
1698 tofid = ( *iter2 )->tofId();
1699 // Scintillator Endcap TOF
1700 if ( !( ( *iter2 )->is_mrpc() ) )
1701 {
1702 if ( tofid < 48 ) { barrelid = 0; }
1703 if ( tofid > 47 ) { barrelid = 2; }
1704 if ( barrelid == 2 ) { tofid = tofid - 48; }
1705 }
1706 // MRPC Endcap TOF
1707 else if ( ( *iter2 )->is_mrpc() )
1708 {
1709 if ( ( TofID::endcap( Identifier( ( *iter2 )->identify() ) ) ) == 0 ) { barrelid = 0; }
1710 else if ( ( TofID::endcap( Identifier( ( *iter2 )->identify() ) ) ) == 1 )
1711 { barrelid = 2; }
1712 if ( barrelid == 2 ) { tofid = tofid - 36; }
1713 }
1714
1715 log << MSG::INFO << " is_mrpc = " << ( *iter2 )->is_mrpc() << " TofId = " << tofid
1716 << " barrelid = " << barrelid << endmsg << " ForwordADC = " << ( *iter2 )->adc()
1717 << " ForwordTDC = " << ( *iter2 )->tdc() << endmsg;
1718 double ftdc = ( *iter2 )->tdc();
1719 double fadc = ( *iter2 )->adc();
1720 if ( m_debug == 4 )
1721 cout << "endcap::single hit,barrelid,tofid,tdc: " << barrelid << " , " << tofid
1722 << " , " << ftdc << endl;
1723
1724 // Scintillator Endcap TOF
1725 if ( !( ( *iter2 )->is_mrpc() ) && useEtofScin )
1726 {
1727 int idptof = ( ( tofid - 1 ) == -1 ) ? 47 : tofid - 1;
1728 int idntof = ( ( tofid + 1 ) == 48 ) ? 0 : tofid + 1;
1729 // Collision Case
1730 if ( idmatch[barrelid][tofid] == 1 || idmatch[barrelid][idptof] == 1 ||
1731 idmatch[barrelid][idntof] == 1 )
1732 {
1733 for ( int i = 0; i < ntot; i++ )
1734 {
1735 if ( ttof[i] != 0 && ftdc > 0 && ftdc < 2000. )
1736 {
1737 if ( ( tofid_helix[i] == tofid ) || ( tofid_helix[i] == idntof ) ||
1738 ( tofid_helix[i] == idptof ) )
1739 {
1740 if ( barrelid == 0 || barrelid == 2 )
1741 {
1742 if ( r_endtof[i] >= 41 && r_endtof[i] <= 90 )
1743 {
1744 if ( optCosmic && ( tofid < 24 || ( tofid > 48 && tofid < 71 ) ) )
1745 {
1746 forevtime = -ttof[i] + r_endtof[i] * 0.09 + 12.2;
1747 meantevup[ntofup] = forevtime;
1748 ntofup++;
1749 }
1750 else
1751 {
1752 forevtime = ttof[i] + r_endtof[i] * 0.09 + 12.2;
1753 meantevdown[ntofdown] = forevtime;
1754 ntofdown++;
1755 }
1756 if ( ( *iter2 )->adc() < 0.0 || m_userawtime )
1757 { t0forward_trk = ftdc - forevtime; }
1758 else
1759 {
1760 t0forward_trk = tofCaliSvc->ETime( ( *iter2 )->adc(), ( *iter2 )->tdc(),
1761 r_endtof[i], ( *iter2 )->tofId() ) -
1762 ttof[i];
1763 }
1764
1765 if ( t0forward_trk < -1. ) continue;
1766 if ( !m_TofOpt && nmatch_end != 0 &&
1767 fabs( t0forward_trk - t0forward_add / nmatch_end ) > 9 )
1768 continue;
1769 t0forward_add += t0forward_trk;
1770 if ( nmatch > 499 ) break;
1771 Tof_t0Arr[nmatch] = t0forward_trk;
1772 meantev[nmatch] = forevtime / 2;
1773 nmatch++;
1774 nmatch_end++;
1775 }
1776 }
1777 }
1778 if ( m_debug == 4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
1779 }
1780 }
1781 }
1782 }
1783 // MRPC Endcap TOF
1784 if ( ( *iter2 )->is_mrpc() && useEtofMRPC )
1785 {
1786 if ( ( ( *iter2 )->quality() & 0x5 ) != 0x5 ) continue;
1787 double btdc = ( *iter2 )->tdc2();
1788 double badc = ( *iter2 )->adc2();
1789 int idptof = ( ( tofid - 1 ) == -1 ) ? 35 : tofid - 1;
1790 int idntof = ( ( tofid + 1 ) == 36 ) ? 0 : tofid + 1;
1791 // B e a m d a t a c a s e <--- Implement for test!
1792 if ( idetfmatch[barrelid][tofid] == 1 || idetfmatch[barrelid][idptof] == 1 ||
1793 idetfmatch[barrelid][idntof] == 1 )
1794 {
1795 for ( int i = 0; i < ntot; i++ )
1796 {
1797 if ( tetf[i] != 0 && ftdc > 0 && ftdc < 2000. )
1798 {
1799 if ( etfid_helix[i] == tofid || etfid_helix[i] == idntof ||
1800 etfid_helix[i] == idptof )
1801 {
1802 if ( barrelid == 0 || barrelid == 2 )
1803 {
1804 if ( r_endetf[i] >= 41 && r_endetf[i] <= 90 )
1805 {
1806 if ( optCosmic && ( tofid < 18 || ( tofid > 35 && tofid < 54 ) ) )
1807 {
1808 forevtime = -tetf[i] + 17.7;
1809 meantevup[ntofup] = forevtime;
1810 ntofup++;
1811 }
1812 else
1813 {
1814 forevtime = tetf[i] + 17.7;
1815 meantevdown[ntofdown] = forevtime;
1816 ntofdown++;
1817 }
1818 if ( m_userawtime )
1819 {
1820 double fbtdc = ( ftdc + btdc ) / 2.0;
1821 if ( ftdc > 0 && btdc < 0 ) { fbtdc = ftdc; }
1822 else if ( ftdc < 0 && btdc > 0 ) { fbtdc = btdc; }
1823 else if ( ftdc < 0 && btdc < 0 ) continue;
1824 t0forward_trk = fbtdc - forevtime;
1825 }
1826 else
1827 {
1828 t0forward_trk =
1829 tofCaliSvc->EtfTime( ( *iter2 )->tdc1(), ( *iter2 )->tdc2(),
1830 ( *iter2 )->tofId(), ( *iter2 )->strip() ) -
1831 tetf[i];
1832 }
1833
1834 if ( t0forward_trk < -1 ) continue;
1835 if ( m_TofOpt && nmatch_end != 0 &&
1836 fabs( t0forward_trk - t0forward_add / nmatch_end ) > 9 )
1837 continue;
1838 if ( m_debug == 4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
1839 t0forward_add += t0forward_trk;
1840 if ( nmatch > 499 ) break;
1841 Tof_t0Arr[nmatch] = t0forward_trk;
1842 meantev[nmatch] = forevtime;
1843 nmatch++;
1844 nmatch_end++;
1845 }
1846 }
1847 }
1848 }
1849 }
1850 }
1851 }
1852 }
1853 if ( nmatch_end ) { tof_flag = 5; }
1854 }
1855
1856 if ( nmatch_barrel == 0 && nmatch_end == 0 )
1857 {
1858 for ( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end();
1859 iter2++, digiId++ )
1860 {
1861 barrelid = ( *iter2 )->barrel();
1862 if ( ( *iter2 )->barrel() != 0 ) continue;
1863 if ( ( *iter2 )->times() != 1 ) continue;
1864 tofid = ( *iter2 )->tofId();
1865 // Scintillator Endcap TOF
1866 if ( !( ( *iter2 )->is_mrpc() ) )
1867 {
1868 if ( tofid < 48 ) { barrelid = 0; }
1869 if ( tofid > 47 ) { barrelid = 2; }
1870 if ( barrelid == 2 ) { tofid = tofid - 48; }
1871 }
1872 // MRPC Endcap TOF
1873 else if ( ( *iter2 )->is_mrpc() )
1874 {
1875 if ( ( TofID::endcap( Identifier( ( *iter2 )->identify() ) ) ) == 0 ) { barrelid = 0; }
1876 else if ( ( TofID::endcap( Identifier( ( *iter2 )->identify() ) ) ) == 1 )
1877 { barrelid = 2; }
1878 if ( barrelid == 2 ) { tofid = tofid - 36; }
1879 }
1880
1881 log << MSG::INFO << " is_mrpc = " << ( *iter2 )->is_mrpc() << " TofId = " << tofid
1882 << " barrelid = " << barrelid << endmsg << " ForwordADC = " << ( *iter2 )->adc()
1883 << " ForwordTDC = " << ( *iter2 )->tdc() << endmsg;
1884 double ftdc = ( *iter2 )->tdc();
1885 double fadc = ( *iter2 )->adc();
1886 if ( m_debug == 4 )
1887 cout << "endcap::single hit,barrelid,tofid,tdc: " << barrelid << " , " << tofid
1888 << " , " << ftdc << endl;
1889
1890 // Scintillator Endcap TOF
1891 if ( !( ( *iter2 )->is_mrpc() ) && useEtofScin )
1892 {
1893 int idptof = ( ( tofid - 1 ) == -1 ) ? 47 : tofid - 1;
1894 int idntof = ( ( tofid + 1 ) == 48 ) ? 0 : tofid + 1;
1895 for ( int i = 0; i < 2; i++ )
1896 {
1897 if ( zemc_rec[0] || zemc_rec[1] )
1898 {
1899 if ( tofid == tofid_emc[i] || tofid_emc[i] == idntof || tofid_emc[i] == idptof )
1900 {
1901 if ( ftdc > 2000. || module[i] == 1 ) continue; // module[i]==1 barrel
1902 ttof[i] =
1903 sqrt( ELMAS2 / ( m_ebeam * m_ebeam ) + 1 ) *
1904 sqrt( xemc_rec[i] * xemc_rec[i] + yemc_rec[i] * yemc_rec[i] + 133 * 133 ) /
1905 VLIGHT;
1906 r_endtof[i] = sqrt( xemc_rec[i] * xemc_rec[i] + yemc_rec[i] * yemc_rec[i] );
1907 if ( optCosmic && ( tofid < 24 || ( tofid > 48 && tofid < 71 ) ) )
1908 {
1909 forevtime = -ttof[i] + r_endtof[i] * 0.09 + 12.2;
1910 meantevup[ntofup] = forevtime;
1911 ntofup++;
1912 }
1913 else
1914 {
1915 forevtime = ttof[i] + r_endtof[i] * 0.09 + 12.2;
1916 meantevdown[ntofdown] = forevtime;
1917 ntofdown++;
1918 }
1919 if ( ( *iter2 )->adc() < 0.0 || m_userawtime )
1920 { t0forward_trk = ftdc - forevtime; }
1921 else
1922 {
1923 t0forward_trk = tofCaliSvc->ETime( ( *iter2 )->adc(), ( *iter2 )->tdc(),
1924 r_endtof[i], ( *iter2 )->tofId() ) -
1925 ttof[i];
1926 if ( m_debug == 4 )
1927 cout << "emc flag t0forward_trk: " << t0forward_trk << endl;
1928 }
1929
1930 if ( t0forward_trk < -1. ) continue;
1931 if ( !m_TofOpt && nmatch_end != 0 &&
1932 fabs( t0forward_trk - t0forward_add / nmatch_end ) > 9 )
1933 continue;
1934 t0forward_add += t0forward_trk;
1935 if ( nmatch > 499 ) break;
1936 meantev[nmatch] = forevtime;
1937 Tof_t0Arr[nmatch] = t0forward_trk;
1938 nmatch++;
1939 nmatch_end++;
1940 emcflag2 = 1;
1941 }
1942 }
1943 }
1944 }
1945 // MRPC Endcap TOF
1946 if ( ( *iter2 )->is_mrpc() && useEtofMRPC )
1947 {
1948 double btdc = ( *iter2 )->tdc2();
1949 double badc = ( *iter2 )->adc2();
1950 int idptof = ( ( tofid - 1 ) == -1 ) ? 35 : tofid - 1;
1951 int idntof = ( ( tofid + 1 ) == 36 ) ? 0 : tofid + 1;
1952 for ( int i = 0; i < 2; i++ )
1953 {
1954 if ( zemc_rec[0] || zemc_rec[1] )
1955 {
1956 if ( tofid == etfid_emc[i] || etfid_emc[i] == idntof || etfid_emc[i] == idptof )
1957 {
1958
1959 if ( ftdc > 2000. || module[i] == 1 ) continue;
1960 tetf[i] =
1961 sqrt( ELMAS2 / ( m_ebeam * m_ebeam ) + 1 ) *
1962 sqrt( xemc_rec[i] * xemc_rec[i] + yemc_rec[i] * yemc_rec[i] + 133 * 133 ) /
1963 VLIGHT;
1964 r_endetf[i] = sqrt( xemc_rec[i] * xemc_rec[i] + yemc_rec[i] * yemc_rec[i] );
1965 if ( optCosmic && ( tofid < 18 || ( tofid > 35 && tofid < 54 ) ) )
1966 {
1967 forevtime = -tetf[i] + 17.7;
1968 meantevup[ntofup] = forevtime;
1969 ntofup++;
1970 }
1971 else
1972 {
1973 forevtime = tetf[i] + 17.7;
1974 meantevdown[ntofdown] = forevtime;
1975 ntofdown++;
1976 }
1977
1978 if ( m_userawtime )
1979 {
1980 double fbtdc = ( ftdc + btdc ) / 2.0;
1981 if ( ftdc > 0 && btdc < 0 ) { fbtdc = ftdc; }
1982 else if ( ftdc < 0 && btdc > 0 ) { fbtdc = btdc; }
1983 else if ( ftdc < 0 && btdc < 0 ) continue;
1984 t0forward_trk = fbtdc - forevtime;
1985 }
1986 else
1987 {
1988 t0forward_trk =
1989 tofCaliSvc->EtfTime( ( *iter2 )->tdc1(), ( *iter2 )->tdc2(),
1990 ( *iter2 )->tofId(), ( *iter2 )->strip() ) -
1991 tetf[i];
1992 }
1993
1994 if ( t0forward_trk < -1 ) continue;
1995 if ( !m_TofOpt && nmatch_end != 0 &&
1996 fabs( t0forward_trk - t0forward_add / nmatch_end ) > 9 )
1997 continue;
1998 if ( m_debug == 4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
1999 t0forward_add += t0forward_trk;
2000 if ( nmatch > 499 ) break;
2001 Tof_t0Arr[nmatch] = t0forward_trk;
2002 nmatch++;
2003 nmatch_end++;
2004 emcflag2 = 1;
2005 }
2006 }
2007 }
2008 }
2009 }
2010 if ( nmatch_end ) { tof_flag = 5; }
2011 }
2012
2013 if ( nmatch_barrel == 0 && nmatch_end == 0 )
2014 {
2015 for ( TofDataVector::iterator iter2 = tofDigiVec.begin(); iter2 != tofDigiVec.end();
2016 iter2++, digiId++ )
2017 {
2018 log << MSG::INFO << "TOF digit No: " << digiId << endmsg;
2019 barrelid = ( *iter2 )->barrel();
2020 if ( ( *iter2 )->barrel() != 0 ) continue;
2021 if ( ( *iter2 )->times() > 1 && ( *iter2 )->times() < 5 )
2022 {
2023 tofid = ( *iter2 )->tofId();
2024 // Scintillator Endcap TOF
2025 if ( !( ( *iter2 )->is_mrpc() ) )
2026 {
2027 if ( tofid < 48 ) { barrelid = 0; }
2028 if ( tofid > 47 ) { barrelid = 2; }
2029 if ( barrelid == 2 ) { tofid = tofid - 48; }
2030 }
2031 // MRPC Endcap TOF
2032 else if ( ( *iter2 )->is_mrpc() )
2033 {
2034 if ( ( TofID::endcap( Identifier( ( *iter2 )->identify() ) ) ) == 0 )
2035 { barrelid = 0; }
2036 else if ( ( TofID::endcap( Identifier( ( *iter2 )->identify() ) ) ) == 1 )
2037 { barrelid = 2; }
2038 if ( barrelid == 2 ) { tofid = tofid - 36; }
2039 }
2040 log << MSG::INFO << " TofId = " << tofid << " barrelid = " << barrelid << endmsg
2041 << " ForwordADC = " << ( *iter2 )->adc() << " ForwordTDC = " << ( *iter2 )->tdc()
2042 << endmsg;
2043 double ftdc = ( *iter2 )->tdc();
2044 double fadc = ( *iter2 )->adc();
2045 if ( m_debug == 4 )
2046 {
2047 cout << "endcap::multi hit,barrelid,tofid,tdc: " << barrelid << " , " << tofid
2048 << " , " << ftdc << endl;
2049 }
2050
2051 // Scintillator Endcap TOF
2052 if ( !( ( *iter2 )->is_mrpc() ) && useEtofScin )
2053 {
2054 int idptof = ( ( tofid - 1 ) == -1 ) ? 47 : tofid - 1;
2055 int idntof = ( ( tofid + 1 ) == 48 ) ? 0 : tofid + 1;
2056 // B e a m d a t a c a s e
2057 if ( idmatch[barrelid][tofid] == 1 || idmatch[barrelid][idptof] == 1 ||
2058 idmatch[barrelid][idntof] == 1 )
2059 {
2060 for ( int i = 0; i < ntot; i++ )
2061 {
2062 if ( ttof[i] != 0 && ftdc > 0 )
2063 {
2064 if ( ( tofid_helix[i] == tofid ) || ( tofid_helix[i] == idntof ) ||
2065 ( tofid_helix[i] == idptof ) )
2066 {
2067 if ( barrelid == 0 || barrelid == 2 )
2068 {
2069 if ( r_endtof[i] >= 41 && r_endtof[i] <= 90 )
2070 {
2071 if ( optCosmic && ( tofid < 24 || ( tofid > 48 && tofid < 71 ) ) )
2072 {
2073 forevtime = -ttof[i] + r_endtof[i] * 0.09 + 12.2;
2074 meantevup[ntofup] = forevtime;
2075 ntofup++;
2076 }
2077 else
2078 {
2079 forevtime = ttof[i] + r_endtof[i] * 0.09 + 12.2;
2080 meantevdown[ntofdown] = forevtime;
2081 ntofdown++;
2082 }
2083 if ( ( *iter2 )->adc() < 0.0 || m_userawtime )
2084 { t0forward_trk = ftdc - forevtime; }
2085 else
2086 {
2087 t0forward_trk =
2088 tofCaliSvc->ETime( ( *iter2 )->adc(), ( *iter2 )->tdc(),
2089 r_endtof[i], ( *iter2 )->tofId() ) -
2090 ttof[i];
2091 }
2092
2093 if ( t0forward_trk < -1. ) continue;
2094 if ( !m_TofOpt && nmatch_end != 0 &&
2095 fabs( t0forward_trk - t0forward_add / nmatch_end ) > 9 )
2096 continue;
2097 t0forward_add += t0forward_trk;
2098 if ( nmatch > 499 ) break;
2099 meantev[nmatch] = forevtime;
2100 Tof_t0Arr[nmatch] = t0forward_trk;
2101 nmatch++;
2102 nmatch_end++;
2103 }
2104 }
2105 if ( m_debug == 4 ) { cout << "t0forward_trk:" << t0forward_trk << endl; }
2106 }
2107 }
2108 }
2109 }
2110 }
2111 // MRPC Endcap TOF
2112 if ( ( *iter2 )->is_mrpc() && useEtofMRPC )
2113 {
2114 double btdc = ( *iter2 )->tdc2();
2115 double badc = ( *iter2 )->adc2();
2116 int idptof = ( ( tofid - 1 ) == -1 ) ? 35 : tofid - 1;
2117 int idntof = ( ( tofid + 1 ) == 36 ) ? 0 : tofid + 1;
2118 // B e a m d a t a c a s e <--- Implement for test!
2119 if ( idetfmatch[barrelid][tofid] == 1 || idetfmatch[barrelid][idptof] == 1 ||
2120 idetfmatch[barrelid][idntof] == 1 )
2121 {
2122 for ( int i = 0; i < ntot; i++ )
2123 {
2124 if ( tetf[i] != 0 && ftdc > 0 && ftdc < 2000. )
2125 {
2126 if ( etfid_helix[i] == tofid || etfid_helix[i] == idntof ||
2127 etfid_helix[i] == idptof )
2128 {
2129 if ( barrelid == 0 || barrelid == 2 )
2130 {
2131 if ( r_endetf[i] >= 41 && r_endetf[i] <= 90 )
2132 {
2133 if ( optCosmic && ( tofid < 18 || ( tofid > 35 && tofid < 54 ) ) )
2134 {
2135 forevtime = -tetf[i] + 17.7;
2136 meantevup[ntofup] = forevtime;
2137 ntofup++;
2138 }
2139 else
2140 {
2141 forevtime = tetf[i] + 17.7;
2142 meantevdown[ntofdown] = forevtime;
2143 ntofdown++;
2144 }
2145 if ( m_userawtime )
2146 {
2147 double fbtdc = ( ftdc + btdc ) / 2.0;
2148 if ( ftdc > 0 && btdc < 0 ) { fbtdc = ftdc; }
2149 else if ( ftdc < 0 && btdc > 0 ) { fbtdc = btdc; }
2150 else if ( ftdc < 0 && btdc < 0 ) continue;
2151 t0forward_trk = fbtdc - forevtime;
2152 }
2153 else
2154 {
2155 t0forward_trk =
2156 tofCaliSvc->EtfTime( ( *iter2 )->tdc1(), ( *iter2 )->tdc2(),
2157 ( *iter2 )->tofId(), ( *iter2 )->strip() ) -
2158 tetf[i];
2159 }
2160
2161 if ( t0forward_trk < -1 ) continue;
2162 if ( !m_TofOpt && nmatch_end != 0 &&
2163 fabs( t0forward_trk - t0forward_add / nmatch_end ) > 9 )
2164 continue;
2165 if ( m_debug == 4 )
2166 { cout << "t0forward_trk:" << t0forward_trk << endl; }
2167 t0forward_add += t0forward_trk;
2168 if ( nmatch > 499 ) break;
2169 Tof_t0Arr[nmatch] = t0forward_trk;
2170 meantev[nmatch] = forevtime;
2171 nmatch++;
2172 nmatch_end++;
2173 }
2174 }
2175 }
2176 }
2177 }
2178 }
2179 }
2180 }
2181 }
2182 if ( nmatch_end ) { tof_flag = 7; }
2183 }
2184
2185 if ( m_ntupleflag && m_tuple2 )
2186 {
2187 g_nmatchbarrel = nmatch_barrel;
2188 g_nmatchbarrel_1 = nmatch_barrel_1;
2189 g_nmatchbarrel_2 = nmatch_barrel_2;
2190 g_nmatchend = nmatch_end;
2191 }
2192
2193 if ( nmatch_end != 0 )
2194 {
2195 t0forward = t0forward_add / nmatch_end;
2196 if ( optCosmic == 0 )
2197 {
2198 if ( m_TofOpt )
2199 {
2200 t_Est = EST_Trimmer( Opt_new( Tof_t0Arr, nmatch, m_TofOpt_Cut ), tOffset_e,
2201 toffset_rawtime_e, offset_dt_end, bunchtime );
2202 }
2203 else
2204 {
2205 t_Est =
2206 EST_Trimmer( t0forward, tOffset_e, toffset_rawtime_e, offset_dt_end, bunchtime );
2207 }
2208 /*
2209 cout << "EsTimeAlg: Determine t_est:" << endl;
2210 cout << " tOffset_b " << tOffset_b << endl;
2211 cout << " toffset_rawtime " <<toffset_rawtime<< endl;
2212 cout << " offset_dt "<< offset_dt << endl;
2213 cout << " Opt_new "<< Opt_new(Tof_t0Arr,nmatch,m_TofOpt_Cut) << endl;
2214 cout << " Tof_t0Arr "<< *Tof_t0Arr << endl;
2215 cout << " nmatch "<< nmatch << endl;
2216 cout << " t_Est " << t_Est << endl;
2217 cout << " t_0forward " << t0forward << endl;
2218 cout << " tOffset_e " << tOffset_e << endl;
2219 cout << " toffset_raw_e " << toffset_rawtime_e << endl;
2220 cout << " offset_dt_end " << offset_dt_end << endl;
2221 cout << " bunchtime " << bunchtime << endl;
2222 cout << " m_TofOpt " << m_TofOpt << endl;
2223 cout << "--------------------------" << endl;
2224 */
2225 if ( t_Est < 0 ) t_Est = 0;
2226 if ( tof_flag == 5 ) tEstFlag = 151;
2227 else if ( tof_flag == 7 ) tEstFlag = 171;
2228 if ( emcflag2 == 1 ) tEstFlag = 161;
2229 /* if(m_nbunch==6){
2230 nbunch=(int)(t0forward-offset)/4;
2231 if(((int)(t0forward-offset))%4>2 || ((int)(t0forward-offset))%4==2) nbunch=nbunch+1;
2232 // if(((int)(t0forward-offset))%4>2) nbunch=nbunch+1;
2233 t_Est=nbunch*4+offset;
2234 if(t_Est<0) t_Est=0;
2235 tEstFlag=151;
2236 }*/
2237 }
2238 if ( optCosmic )
2239 {
2240 t_Est = t0forward; // 10. yzhang add for nmatch_end cosmic event
2241 if ( tof_flag == 5 ) tEstFlag = 251;
2242 else if ( tof_flag == 7 ) tEstFlag = 271;
2243 if ( emcflag2 == 1 ) tEstFlag = 261;
2244 }
2245 if ( m_ntupleflag && m_tuple2 ) g_meant0 = t0forward;
2246 }
2247
2248 double t0_comp = -999;
2249 double T0 = -999;
2250
2251 if ( nmatch_barrel == 0 && nmatch_end == 0 && m_flag == 1 )
2252 {
2253 double mhit[43][300] = { 0. };
2254 SmartDataPtr<MdcDigiCol> mdcDigiCol( eventSvc(), "/Event/Digi/MdcDigiCol" );
2255 if ( !mdcDigiCol )
2256 {
2257 log << MSG::INFO << "Could not find MDC digi" << endmsg;
2258 return StatusCode::FAILURE;
2259 }
2260
2261 IMdcGeomSvc* mdcGeomSvc;
2262 StatusCode sc = service( "MdcGeomSvc", mdcGeomSvc );
2263 if ( sc != StatusCode::SUCCESS ) { return StatusCode::FAILURE; }
2264
2265 MdcDigiCol::iterator iter1 = mdcDigiCol->begin();
2266 digiId = 0;
2267 Identifier mdcId;
2268 int layerId;
2269 int wireId;
2270
2271 for ( ; iter1 != mdcDigiCol->end(); iter1++, digiId++ )
2272 {
2273 mdcId = ( *iter1 )->identify();
2274 layerId = MdcID::layer( mdcId );
2275 wireId = MdcID::wire( mdcId );
2276
2277 mhit[layerId][wireId] = RawDataUtil::MdcTime( ( *iter1 )->getTimeChannel() );
2278 mhit[layerId][wireId] -= 1.28 * ( mdcGeomSvc->Layer( layerId )->Radius() ) / 299.8;
2279
2280 mdcGeomSvc->Wire( layerId, wireId );
2281 // Apply crude TOF, Belle: tof=1.074*radius/c, here we use 1.28 instead of 1.074.
2282 double tof;
2283 // tof = 1.28 * mhit.geo->Lyr()->Radius()/299.8; //unit of Radius is mm.
2284 }
2285
2286 int Iused[43][300] = { 0 }, tmp = 0;
2287 bool Lpat, Lpat11, Lpat12, Lpat2, Lpat31, Lpat32;
2288 double t0_all = 0, t0_mean = 0;
2289 double r[4] = { 0. };
2290 double chi2 = 999.0;
2291 double phi[4] = { 0. }, corr[4] = { 0. }, driftt[4] = { 0. };
2292 double ambig = 1;
2293 double mchisq = 50000;
2294
2295 // for(int i=2;i<10;i++)
2296 for ( int i = 5; i < 10; i++ )
2297 {
2298
2299 double T1 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * i + 0 )->PCSiz() ) / 0.004;
2300 double T2 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * i + 1 )->PCSiz() ) / 0.004;
2301 double T3 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * i + 2 )->PCSiz() ) / 0.004;
2302 double T4 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * i + 3 )->PCSiz() ) / 0.004;
2303 r[0] = ( mdcGeomSvc->Layer( 4 * i + 0 )->Radius() ) * 0.1;
2304 r[1] = ( mdcGeomSvc->Layer( 4 * i + 1 )->Radius() ) * 0.1;
2305 r[2] = ( mdcGeomSvc->Layer( 4 * i + 2 )->Radius() ) * 0.1;
2306 r[3] = ( mdcGeomSvc->Layer( 4 * i + 3 )->Radius() ) * 0.1;
2307 double r0 = r[0] - r[1] - ( r[2] - r[1] ) / 2;
2308 double r1 = -( r[2] - r[1] ) / 2;
2309 double r2 = ( r[2] - r[1] ) / 2;
2310 double r3 = r[3] - r[2] + ( r[2] - r[1] ) / 2;
2311
2312 for ( int j = 0; j < mdcGeomSvc->Layer( i * 4 + 3 )->NCell(); j++ )
2313 {
2314
2315 int Icp = 0;
2316 Icp = j - 1;
2317 if ( Icp < 0 ) Icp = mdcGeomSvc->Layer( i * 4 + 3 )->NCell();
2318
2319 Lpat = ( mhit[4 * i][j] != 0 ) && ( mhit[4 * i][Icp] == 0 ) &&
2320 ( mhit[4 * i][j + 1] == 0 ) && ( Iused[4 * i][j] == 0 );
2321 if ( Lpat )
2322 {
2323 Lpat11 = ( mhit[4 * i + 1][Icp] == 0 ) && ( Iused[4 * i + 1][j] == 0 ) &&
2324 ( mhit[4 * i + 1][j] != 0 ) && ( mhit[4 * i + 1][j + 1] == 0 );
2325 Lpat12 = ( mhit[4 * i + 1][j] == 0 ) && ( Iused[4 * i + 1][j + 1] == 0 ) &&
2326 ( mhit[4 * i + 1][j + 1] != 0 ) && ( mhit[4 * i + 1][j + 2] == 0 );
2327 Lpat2 = ( mhit[4 * i + 2][Icp] == 0 ) && ( Iused[4 * i + 2][j] == 0 ) &&
2328 ( mhit[4 * i + 2][j] != 0 ) && ( mhit[4 * i + 2][j + 1] == 0 );
2329 Lpat31 = ( mhit[4 * i + 3][Icp] == 0 ) && ( Iused[4 * i + 3][j] == 0 ) &&
2330 ( mhit[4 * i + 3][j] != 0 ) && ( mhit[4 * i + 3][j + 1] == 0 );
2331 Lpat32 = ( mhit[4 * i + 3][j] == 0 ) && ( Iused[4 * i + 3][j + 1] == 0 ) &&
2332 ( mhit[4 * i + 3][j + 1] != 0 ) && ( mhit[4 * i + 3][j + 2] == 0 );
2333
2334 if ( Lpat11 && Lpat2 && Lpat31 )
2335 {
2336
2337 Iused[4 * i + 0][j] = 1;
2338 Iused[4 * i + 1][j] = 1;
2339 Iused[4 * i + 2][j] = 1;
2340 Iused[4 * i + 3][j] = 1;
2341 double t_i = mhit[4 * i + 0][j] + mhit[4 * i + 2][j];
2342 double t_j = mhit[4 * i + 1][j] + mhit[4 * i + 3][j];
2343 double l_j = T2 + T4;
2344 double r_i = r0 + r2;
2345 double r_j = r1 + r3;
2346 double r_2k = r0 * r0 + r1 * r1 + r2 * r2 + r3 * r3;
2347 double rt_i = r0 * mhit[4 * i + 0][j] + r2 * mhit[4 * i + 2][j];
2348 double rt_j = r1 * mhit[4 * i + 1][j] + r3 * mhit[4 * i + 3][j];
2349 double rl_j = r1 * T2 + r3 * T4;
2350
2351 double deno = 4 * r_2k - 2 * ( r_i * r_i + r_j * r_j );
2352
2353 if ( deno != 0 )
2354 {
2355 double Pa = ( 4 * ( rt_i - rt_j + rl_j ) - ( t_i + t_j - l_j ) * ( r_i - r_j ) -
2356 ( t_i - t_j + l_j ) * ( r_i + r_j ) ) /
2357 deno;
2358 double Pb = 0.25 * ( t_i - t_j + l_j - ( r_i + r_j ) * Pa );
2359 double Ang = fabs( 90. - fabs( atan( Pa ) * 180. / 3.141593 ) );
2360
2361 t0_all += ( -0.25 * ( ( r_i - r_j ) * Pa - t_i - t_j + l_j ) );
2362
2363 double chi2_tmp;
2364 for ( int t0c = 0; t0c < 17; t0c += 8 )
2365 {
2366 chi2_tmp = ( mhit[4 * i + 0][j] - t0c - r0 * Pa - Pb ) *
2367 ( mhit[4 * i + 0][j] - t0c - r0 * Pa - Pb ) +
2368 ( T2 - mhit[4 * i + 1][j] + t0c - r1 * Pa - Pb ) *
2369 ( T2 - mhit[4 * i + 1][j] + t0c - r1 * Pa - Pb ) +
2370 ( mhit[4 * i + 2][j] - t0c - r2 * Pa - Pb ) *
2371 ( mhit[4 * i + 2][j] - t0c - r2 * Pa - Pb ) +
2372 ( T4 - mhit[4 * i + 3][j] + t0c - r3 * Pa - Pb ) *
2373 ( T4 - mhit[4 * i + 3][j] + t0c - r3 * Pa - Pb );
2374 if ( chi2_tmp < chi2 )
2375 {
2376 chi2 = chi2_tmp;
2377 t0_comp = t0c;
2378 }
2379 }
2380 tmp += 1;
2381 }
2382 // use squareleas t0
2383
2384 for ( int tmpT0 = 0; tmpT0 < 17; tmpT0 += 8 )
2385 {
2386 driftt[0] = mhit[4 * i + 0][j] - tmpT0;
2387 driftt[1] = mhit[4 * i + 1][j] - tmpT0;
2388 driftt[2] = mhit[4 * i + 2][j] - tmpT0;
2389 driftt[3] = mhit[4 * i + 3][j] - tmpT0;
2390
2391 phi[0] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 )->NCell() ) +
2392 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 1 )->NCell() ) / 2;
2393 phi[1] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 1 )->NCell() );
2394 phi[2] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 2 )->NCell() ) +
2395 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 1 )->NCell() ) / 2;
2396 phi[3] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 3 )->NCell() );
2397 phi[0] -= ambig * driftt[0] * 0.004 / r[0];
2398 phi[1] += ambig * driftt[1] * 0.004 / r[1];
2399 phi[2] -= ambig * driftt[2] * 0.004 / r[2];
2400 phi[3] += ambig * driftt[3] * 0.004 / r[3];
2401 double s1, sx, sy, sxx, sxy; // The various sums.
2402 double delinv, denom;
2403 double weight; // weight for hits, calculated from sigmas.
2404 double sigma;
2405 double x[4];
2406 x[0] = r0;
2407 x[1] = r1;
2408 x[2] = r2;
2409 x[3] = r3;
2410 sigma = 0.12;
2411 s1 = sx = sy = sxx = sxy = 0.0;
2412 double chisq = 0.0;
2413 for ( int ihit = 0; ihit < 4; ihit++ )
2414 {
2415 weight = 1. / ( sigma * sigma ); // NEED sigma of MdcHits
2416 s1 += weight;
2417 sx += x[ihit] * weight;
2418 sy += phi[ihit] * weight;
2419 sxx += x[ihit] * ( x[ihit] * weight );
2420 sxy += phi[ihit] * ( x[ihit] * weight );
2421 }
2422 double resid[4] = { 0. };
2423 /* Calculate parameters. */
2424 denom = s1 * sxx - sx * sx;
2425 delinv = ( denom == 0.0 ) ? 1.e20 : 1. / denom;
2426 double intercept = ( sy * sxx - sx * sxy ) * delinv;
2427 double slope = ( s1 * sxy - sx * sy ) * delinv;
2428
2429 /* Calculate residuals. */
2430 for ( int ihit = 0; ihit < 4; ihit++ )
2431 {
2432 resid[ihit] = ( phi[ihit] - intercept - slope * x[ihit] );
2433 chisq += resid[ihit] * resid[ihit] / ( sigma * sigma );
2434 }
2435 if ( chisq < mchisq )
2436 {
2437 mchisq = chisq;
2438 T0 = tmpT0;
2439 }
2440 }
2441 }
2442 if ( Lpat12 && Lpat2 && Lpat32 )
2443 {
2444 Iused[4 * i + 0][j] = 1;
2445 Iused[4 * i + 1][j + 1] = 1;
2446 Iused[4 * i + 2][j] = 1;
2447 Iused[4 * i + 3][j + 1] = 1;
2448
2449 double t_i = mhit[4 * i + 0][j] + mhit[4 * i + 2][j];
2450 double t_j = mhit[4 * i + 1][j + 1] + mhit[4 * i + 3][j + 1];
2451 double l_j = T2 + T4;
2452 double r_i = r0 + r2;
2453 double r_j = r1 + r3;
2454 double r_2k = r0 * r0 + r1 * r1 + r2 * r2 + r3 * r3;
2455 double rt_i = r0 * mhit[4 * i + 0][j] + r2 * mhit[4 * i + 2][j];
2456 double rt_j = r1 * mhit[4 * i + 1][j + 1] + r3 * mhit[4 * i + 3][j + 1];
2457 double rl_j = r1 * T2 + r3 * T4;
2458 double deno = 4 * r_2k - 2 * ( r_i * r_i + r_j * r_j );
2459
2460 if ( deno != 0 )
2461 {
2462 double Pa = ( 4 * ( rt_i - rt_j + rl_j ) - ( t_i + t_j - l_j ) * ( r_i - r_j ) -
2463 ( t_i - t_j + l_j ) * ( r_i + r_j ) ) /
2464 deno;
2465 double Pb = 0.25 * ( t_i - t_j + l_j - ( r_i + r_j ) * Pa );
2466 double Ang = fabs( 90. - fabs( atan( Pa ) * 180. / 3.141593 ) );
2467 t0_all += ( -0.25 * ( ( r_i - r_j ) * Pa - t_i - t_j + l_j ) );
2468 tmp += 1;
2469 double chi2_tmp;
2470
2471 for ( int t0c = 0; t0c < 17; t0c += 8 )
2472 {
2473 chi2_tmp = ( mhit[4 * i + 0][j] - t0c - r0 * Pa - Pb ) *
2474 ( mhit[4 * i + 0][j] - t0c - r0 * Pa - Pb ) +
2475 ( T2 - mhit[4 * i + 1][j + 1] + t0c - r1 * Pa - Pb ) *
2476 ( T2 - mhit[4 * i + 1][j + 1] + t0c - r1 * Pa - Pb ) +
2477 ( mhit[4 * i + 2][j] - t0c - r2 * Pa - Pb ) *
2478 ( mhit[4 * i + 2][j] - t0c - r2 * Pa - Pb ) +
2479 ( T4 - mhit[4 * i + 3][j + 1] + t0c - r3 * Pa - Pb ) *
2480 ( T4 - mhit[4 * i + 3][j + 1] + t0c - r3 * Pa - Pb );
2481
2482 if ( chi2_tmp < chi2 )
2483 {
2484 chi2 = chi2_tmp;
2485 t0_comp = t0c;
2486 }
2487 }
2488 }
2489
2490 // use squareleast
2491
2492 for ( int tmpT0 = 0; tmpT0 < 17; tmpT0 += 8 )
2493 {
2494 driftt[0] = mhit[4 * i + 0][j] - tmpT0;
2495 driftt[1] = mhit[4 * i + 1][j + 1] - tmpT0;
2496 driftt[2] = mhit[4 * i + 2][j] - tmpT0;
2497 driftt[3] = mhit[4 * i + 3][j + 1] - tmpT0;
2498
2499 phi[0] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 )->NCell() ) +
2500 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 1 )->NCell() ) / 2;
2501 phi[1] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 1 )->NCell() );
2502 phi[2] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 2 )->NCell() ) +
2503 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 1 )->NCell() ) / 2;
2504 phi[3] = j * 2 * 3.14159265 / ( mdcGeomSvc->Layer( i * 4 + 3 )->NCell() );
2505 phi[0] -= ambig * driftt[0] * 0.004 / r[0];
2506 phi[1] += ambig * driftt[1] * 0.004 / r[1];
2507 phi[2] -= ambig * driftt[2] * 0.004 / r[2];
2508 phi[3] += ambig * driftt[3] * 0.004 / r[3];
2509 double s1, sx, sy, sxx, sxy; // The various sums.
2510 double delinv, denom;
2511 double weight; // weight for hits, calculated from sigmas.
2512 double sigma;
2513 double x[4];
2514 x[0] = r0;
2515 x[1] = r1;
2516 x[2] = r2;
2517 x[3] = r3;
2518 sigma = 0.12;
2519 s1 = sx = sy = sxx = sxy = 0.0;
2520 double chisq = 0.0;
2521 for ( int ihit = 0; ihit < 4; ihit++ )
2522 {
2523 weight = 1. / ( sigma * sigma ); // NEED sigma of MdcHits
2524 s1 += weight;
2525 sx += x[ihit] * weight;
2526 sy += phi[ihit] * weight;
2527 sxx += x[ihit] * ( x[ihit] * weight );
2528 sxy += phi[ihit] * ( x[ihit] * weight );
2529 }
2530 double resid[4] = { 0. };
2531 /* Calculate parameters. */
2532 denom = s1 * sxx - sx * sx;
2533 delinv = ( denom == 0.0 ) ? 1.e20 : 1. / denom;
2534 double intercept = ( sy * sxx - sx * sxy ) * delinv;
2535 double slope = ( s1 * sxy - sx * sy ) * delinv;
2536
2537 /* Calculate residuals. */
2538 for ( int ihit = 0; ihit < 4; ihit++ )
2539 {
2540 resid[ihit] = ( phi[ihit] - intercept - slope * x[ihit] );
2541 chisq += resid[ihit] * resid[ihit] / ( sigma * sigma );
2542 }
2543
2544 if ( chisq < mchisq )
2545 {
2546 mchisq = chisq;
2547 T0 = tmpT0;
2548 }
2549 }
2550 }
2551 }
2552 }
2553 }
2554
2555 if ( tmp != 0 ) { t0_mean = t0_all / tmp; }
2556 if ( m_ntupleflag && m_tuple2 ) g_t0mean = t0_mean;
2557
2558 t_Est = T0 + tOffset_b; // 11.yzhang add tOffset_b, MDC method calc. Tes
2559 tEstFlag = 2;
2560 }
2561 if ( m_ntupleflag && m_tuple2 )
2562 {
2563 g_t0 = t0_comp;
2564 g_T0 = T0;
2565 }
2566 if ( nmatch_barrel == 0 && nmatch_end == 0 && nmatch_barrel_1 == 0 && nmatch_barrel_2 == 0 &&
2567 m_mdcCalibFunSvc && m_flag == 2 )
2568 {
2569
2570 log << MSG::INFO << " mdc " << endmsg;
2571
2572#define MXWIRE 6860
2573#define MXTKHIT 6860
2574#define MXTRK 15
2575
2576 // L o c a l v a r i a b l e s
2577 int nhits_ax = 0;
2578 int nhits_ax2 = 0;
2579 int nhits_st = 0;
2580 int nhits_st2 = 0;
2581
2582 double tev_ax[MXTKHIT];
2583 double tev_st[MXTKHIT];
2584 double tevv[MXTKHIT];
2585 double toft;
2586 double prop;
2587 double t0_minus_TDC[MXWIRE];
2588 // double adc[MXWIRE];
2589 double T0_cal = -230;
2590 double Mdc_t0Arr[500];
2591
2592 int expmc = 1;
2593 double scale = 1.;
2594 int expno, runno;
2595 ndriftt = 0;
2596
2597 // A l l t r a c k s
2598 // Check if Fzisan track exist
2599 // if(ntot<=0) return 0; // it was "if(ntot<=0) return 0;"
2600 if ( ntot > MXTRK ) { ntot = MXTRK; }
2601
2602 // SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
2603 if ( !newtrkCol || newtrkCol->size() == 0 )
2604 {
2605 log << MSG::INFO << "Could not find RecMdcTrackCol" << endmsg;
2606 return StatusCode::SUCCESS;
2607 }
2608 log << MSG::INFO << "Begin to check RecMdcTrackCol" << endmsg;
2609
2610 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
2611
2612 for ( int tempntrack = 0; iter_trk != newtrkCol->end(); iter_trk++, tempntrack++ )
2613 {
2614 log << MSG::DEBUG << "retrieved MDC track:"
2615 << " Track Id: " << ( *iter_trk )->trackId() << " Dr: " << ( *iter_trk )->helix( 0 )
2616 << " Phi0: " << ( *iter_trk )->helix( 1 ) << " kappa: " << ( *iter_trk )->helix( 2 )
2617 << " Dz: " << ( *iter_trk )->helix( 3 ) << " Tanl: " << ( *iter_trk )->helix( 4 )
2618 << " Phi terminal: " << ( *iter_trk )->getFiTerm() << endmsg
2619 << "Number of hits: " << ( *iter_trk )->getNhits() << " Number of stereo hits "
2620 << ( *iter_trk )->nster() << endmsg;
2621
2622 // Track variables
2623 const HepPoint3D pivot0( 0., 0., 0. );
2624 HepVector a( 5, 0 );
2625
2626 a[0] = ( *iter_trk )->helix( 0 );
2627 a[1] = ( *iter_trk )->helix( 1 );
2628 a[2] = ( *iter_trk )->helix( 2 );
2629 a[3] = ( *iter_trk )->helix( 3 );
2630 a[4] = ( *iter_trk )->helix( 4 );
2631
2632 // Ill-fitted (dz=tanl=-9999.) or off-IP track in fzisan
2633 if ( abs( a[0] ) > Estparam.MDC_drCut() || abs( a[3] ) > Estparam.MDC_dzCut() ||
2634 abs( a[4] ) > 500. )
2635 continue;
2636
2637 double phi0 = a[1];
2638 double kappa = abs( a[2] );
2639 double dirmag = sqrt( 1. + a[4] * a[4] );
2640 // double unit_s = abs(rho * dirmag);
2641 double mom = abs( dirmag / kappa );
2642 double beta = mom / sqrt( mom * mom + PIMAS2 );
2643 if ( particleId[tempntrack] == 1 ) { beta = mom / sqrt( mom * mom + ELMAS2 ); }
2644 if ( particleId[tempntrack] == 5 ) { beta = mom / sqrt( mom * mom + PROTONMAS2 ); }
2645
2646 // D e f i n e h e l i x
2647 Helix helix0( pivot0, a );
2648 double rho = helix0.radius();
2649 double unit_s = abs( rho * dirmag );
2650
2651 helix0.ignoreErrorMatrix();
2652 HepPoint3D hcen = helix0.center();
2653 double xc = hcen( 0 );
2654 double yc = hcen( 1 );
2655
2656 if ( xc == 0.0 && yc == 0.0 ) continue;
2657
2658 double direction = 1.;
2659 if ( optCosmic != 0 )
2660 {
2661 double phi = atan2( helix0.momentum( 0. ).y(), helix0.momentum( 0. ).x() );
2662 if ( phi > 0. && phi <= M_PI ) direction = -1.;
2663 }
2664
2665 IMdcGeomSvc* mdcGeomSvc;
2666 StatusCode sc = service( "MdcGeomSvc", mdcGeomSvc );
2667 double zeast;
2668 double zwest;
2669 double m_vp[43] = { 0. }, m_zst[43] = { 0. };
2670 for ( int lay = 0; lay < 43; lay++ )
2671 {
2672 zwest = mdcGeomSvc->Wire( lay, 0 )->Forward().z();
2673 zeast = mdcGeomSvc->Wire( lay, 0 )->Backward().z();
2674 // m_zwid[lay] = (zeast - zwest) / (double)m_nzbin;
2675
2676 if ( lay < 8 ) m_vp[lay] = 220.0; // *10^9 mm/s
2677 else m_vp[lay] = 240.0; // *10^9 mm/s
2678
2679 if ( 0 == ( lay % 2 ) )
2680 { // west end
2681 m_zst[lay] = zwest;
2682 }
2683 else
2684 { // east end
2685 m_zst[lay] = zeast;
2686 }
2687 }
2688
2689 // Hits
2690 log << MSG::DEBUG << "hitList of this track:" << endmsg;
2691 HitRefVec gothits = ( *iter_trk )->getVecHits();
2692 HitRefVec::iterator it_gothit = gothits.begin();
2693 for ( ; it_gothit != gothits.end(); it_gothit++ )
2694 {
2695
2696 log << MSG::DEBUG << "hits Id: " << ( *it_gothit )->getId()
2697 << " hits MDC layerId wireId " << MdcID::layer( ( *it_gothit )->getMdcId() )
2698 << " " << MdcID::wire( ( *it_gothit )->getMdcId() ) << endmsg << " hits TDC "
2699 << ( *it_gothit )->getTdc() << endmsg;
2700
2701 int layer = MdcID::layer( ( *it_gothit )->getMdcId() );
2702 int wid = MdcID::wire( ( *it_gothit )->getMdcId() );
2703 double tdc = ( *it_gothit )->getTdc();
2704 // cout<<"-----------mdc layer,wireid,tdc--------------:
2705 // "<<layer<<","<<wid<<","<<tdc<<endl;
2706 double trkchi2 = ( *iter_trk )->chi2();
2707 if ( trkchi2 > 100 ) continue;
2708 double hitChi2 = ( *it_gothit )->getChisqAdd();
2709 HepVector helix_par = ( *iter_trk )->helix();
2710 HepSymMatrix helixErr = ( *iter_trk )->err();
2711 // *** A x i a l s e g m e n t s
2712 if ( ( layer >= 8 && layer <= 19 ) || ( layer >= 36 && layer <= 42 ) )
2713 {
2714 // H i t s
2715 ////Geomdc_wire &GeoRef = GeoMgr[wid-1];
2716 // MdcGeoWire & GeoRef = Wire[wid-1];
2717
2718 const MdcGeoWire* GeoRef = mdcGeomSvc->Wire( layer, wid );
2719
2720 ////int layer = GeoRef->Layer();
2721 //// int local = GeoRef->Cell();
2722
2723 // int layer = mdcGeomSvc->Wire(wid-1)->Layer();
2724 // int local = mdcGeomSvc->Wire(wid-1)->Cell();
2725 // double fadc = adc[wid];
2726 //// if(mdc_adc_cut_(&expmc, &expno, &runno, &fadc, &layer, &local)==0) continue;
2727
2728 // Use or not to use inner 4 layers (layers with high bg.)
2729 if ( Estparam.MDC_Inner() == 0 && layer <= 3 ) continue;
2730
2731 double xw = GeoRef->Forward().x() / 10; // wire x position at z=0
2732 double yw = GeoRef->Forward().y() / 10; // wire y position at z=0
2733 // move pivot to the wire (slant angle ignored)
2734 HepPoint3D pivot1( xw, yw, 0. );
2735 helix0.pivot( pivot1 );
2736 double zw = helix0.a()[3]; // z position
2737
2738 // T o F c o r r e c t i o n
2739 double dphi = ( -xc * ( xw - xc ) - yc * ( yw - yc ) ) /
2740 sqrt( ( xc * xc + yc * yc ) *
2741 ( ( xw - xc ) * ( xw - xc ) + ( yw - yc ) * ( yw - yc ) ) );
2742 dphi = acos( dphi );
2743 double pathtof = abs( unit_s * dphi );
2744 if ( kappa != 0 ) { toft = pathtof / VLIGHT / beta; }
2745 else { toft = pathtof / VLIGHT; }
2746
2747 // P r o p a g a t i o n d e l a y c o r r e c t i o n
2748 ////if (zw < GeoRef.zwb()) zw = GeoRef.zwb();
2749 ////if (zw > GeoRef.zwf()) zw = GeoRef.zwf();
2750
2751 if ( zw > ( GeoRef->Backward().z() ) / 10 ) zw = ( GeoRef->Backward().z() ) / 10;
2752 if ( zw < ( GeoRef->Forward().z() ) / 10 ) zw = ( GeoRef->Forward().z() ) / 10;
2753
2754 double slant = GeoRef->Slant();
2755 prop = zw / cos( slant ) / VELPROP;
2756 // Time walk correction
2757 double Tw = 0;
2758 // if(expmc==1) {
2759 // Calmdc_const3 &TwRef = Cal3Mgr[layer];
2760 // if(adc[wid]>0.) Tw = TwRef.tw(0) + TwRef.tw(1)/sqrt(adc[wid]);
2761 // }
2762
2763 double driftt;
2764 double dummy;
2765 int lr = 2;
2766 // if((xw*helix0.x(0.).y()-yw*helix0.x(0.).x())<0) lr=-1;
2767 double p[3];
2768 p[0] = helix0.momentum( 0. ).x();
2769 p[1] = helix0.momentum( 0. ).y();
2770 double pos[2];
2771 pos[0] = xw;
2772 pos[1] = yw;
2773 double alpha;
2774 //// calmdc_getalpha_( p, pos, &alpha);
2775 // double dist = wir.distance(); this is wrong
2776 // double dist = abs(helix0.a()[0]); this if wrong
2777 double dist = 0.;
2778
2779 dist = ( m_mdcUtilitySvc->doca( layer, wid, helix_par, helixErr ) ) * 10.0;
2780
2781 if ( dist < 0. ) lr = 1;
2782 else lr = 0;
2783 dist = fabs( dist );
2784 if ( dist > 0.4 * ( mdcGeomSvc->Layer( layer ) )->PCSiz() ) continue;
2785
2786 //* drift distance cut
2787 // int ia;
2788 // ia = (int) ((alpha+90.)/10.);
2789 // double da = alpha+90. - ia*10.;
2790 // if (ia==0) ia=18;
2791
2792 // Calmdc_const &BndRef = Cal1Mgr[layer];
2793 // if(lr==-1) {
2794 // if(dist < BndRef.bndl(ia,0)) continue;
2795 // if(dist > BndRef.bndl(ia,3)) continue;
2796 //}
2797 // if(lr== 1) {
2798 // if(dist < BndRef.bndr(ia,0)) continue;
2799 // if(dist > BndRef.bndr(ia,3)) continue;
2800 // }
2801
2802 int idummy;
2803
2804 if ( !m_useXT ) { driftt = dist / Estparam.vdrift(); }
2805 else
2806 {
2807 double entrance = ( *it_gothit )->getEntra();
2808 driftt = m_mdcCalibFunSvc->distToDriftTime( dist, layer, wid, lr, entrance );
2809 }
2810 if ( m_useT0cal )
2811 {
2812 T0_cal = m_mdcCalibFunSvc->getT0( layer, wid ) +
2813 m_mdcCalibFunSvc->getTimeWalk( layer, tdc );
2814 }
2815
2816 double zprop = fabs( zw - m_zst[layer] );
2817 double tp = zprop / m_vp[layer];
2818 // cout<<"proptation time --tp ax: "<<tp<<endl;
2819 if ( driftt > tdc ) continue;
2820 double difft = tdc - driftt - toft - tp - T0_cal;
2821 if ( ndriftt >= 500 ) break;
2822 if ( difft < -10 ) continue;
2823 Mdc_t0Arr[ndriftt] = difft;
2824 // cout<<"ax: tdc, driftt, toft, tp: "<<tdc<<" , "<<driftt<<" , "<<toft<<" , "<<tp<<"
2825 // , "<<difft<<endl;
2826 sum_EstimeMdc = sum_EstimeMdc + difft;
2827 ndriftt++;
2828 /*if(Estparam.MDC_Xt()==2) {
2829 calmdc_d2t_bfld_( &driftt, &dist, &dummy, &dummy, &lr,
2830 &idummy, &layer, &alpha, &dummy, &dummy, &dummy);
2831 } else if(Estparam.MDC_Xt()==1) {
2832 driftt = dist/Estparam.vdrift();
2833
2834 }*/
2835 // htf[1]->accumulate( t0_minus_TDC[wid] );
2836
2837 double tev = -t0_minus_TDC[wid] + driftt;
2838 if ( Estparam.MDC_Tof() != 0 ) tev += direction * toft;
2839 if ( Estparam.MDC_Prop() != 0 ) tev += prop;
2840 //// if(Estparam.MDC_Walk()!=0) tev += Tw;
2841
2842 // sum_tev_ax += tev;
2843 // htf[3+hid]->accumulate( tev );
2844 nhits_ax++;
2845 tev_ax[nhits_ax - 1] = tev;
2846
2847 if ( Estparam.MDC_Debug() != 0 ) log << MSG::INFO << "*** tev ***" << tev << endmsg;
2848 double driftt_mea = t0_minus_TDC[wid];
2849 // if(abs(driftt-t0_minus_TDC[wid])<75.)
2850 if ( abs( driftt - driftt_mea ) < 75. )
2851 {
2852 // sum_tev_ax2 += tev;
2853 nhits_ax2++;
2854 if ( Estparam.MDC_Debug() != 0 )
2855 log << MSG::INFO << "*** tev2 ***" << tev << endmsg;
2856 }
2857 } // End axial hits
2858
2859 // S t e r e o s e g m e n t s
2860 else if ( ( ( layer >= 4 && layer <= 7 ) || ( layer >= 20 && layer <= 35 ) ) &&
2861 m_useSw )
2862 {
2863
2864 IMdcGeomSvc* mdcGeomSvc;
2865 StatusCode sc = service( "MdcGeomSvc", mdcGeomSvc );
2866 const MdcGeoWire* GeoRef = mdcGeomSvc->Wire( layer, wid );
2867
2868 // double fadc=adc[wid];
2869 //// if(mdc_adc_cut_(&expmc, &expno, &runno, &fadc, &layer, &local)==0) continue;
2870
2871 double bx = GeoRef->Backward().x() / 10;
2872 double by = GeoRef->Backward().y() / 10;
2873 double bz = GeoRef->Backward().z() / 10;
2874 double fx = GeoRef->Forward().x() / 10;
2875 double fy = GeoRef->Forward().y() / 10;
2876 double fz = GeoRef->Forward().z() / 10;
2877
2878 //// HepPoint3D
2879 /// fwd(GeoRef->Forward().x(),GeoRef->Forward().y(),GeoRef->Forward().z()); /
2880 /// HepPoint3D
2881 /// bck(GeoRef->Backward().x(),GeoRef->Backward().y(),GeoRef->Backward().z());
2882 HepPoint3D fwd( fx, fy, fz );
2883 HepPoint3D bck( bx, by, bz );
2884
2885 Hep3Vector wire = (CLHEP::Hep3Vector)bck - (CLHEP::Hep3Vector)fwd;
2886 HepPoint3D try1 = ( fwd + bck ) * .5;
2887 helix0.pivot( try1 );
2888 HepPoint3D try2 = ( helix0.x( 0 ).z() - bck.z() ) / wire.z() * wire + bck;
2889 helix0.pivot( try2 );
2890 HepPoint3D try3 = ( helix0.x( 0 ).z() - bck.z() ) / wire.z() * wire + bck;
2891 helix0.pivot( try3 );
2892
2893 double xh = helix0.x( 0. ).x();
2894 double yh = helix0.x( 0. ).y();
2895 double z = helix0.x( 0. ).z();
2896
2897 // T o F c o r r e c t i o n
2898 double dphi = ( -xc * ( xh - xc ) - yc * ( yh - yc ) ) /
2899 sqrt( ( xc * xc + yc * yc ) *
2900 ( ( xh - xc ) * ( xh - xc ) + ( yh - yc ) * ( yh - yc ) ) );
2901 dphi = acos( dphi );
2902 double pathtof = abs( unit_s * dphi );
2903 if ( kappa != 0 ) { toft = pathtof / VLIGHT / beta; }
2904 else { toft = pathtof / VLIGHT; }
2905
2906 // P r o p a g a t i o n d e l a y c o r r e c t i o n
2907
2908 if ( z != 9999. )
2909 {
2910 // if (z < GeoRef->Forward().z()/10) z = GeoRef->Forward().z()/10;
2911 if ( z < fz ) z = fz;
2912 // if (z >GeoRef->Backward().z()/10) z =GeoRef->Backward().z()/10;
2913 if ( z > bz ) z = bz;
2914 double slant = GeoRef->Slant();
2915 prop = z / cos( slant ) / VELPROP;
2916 }
2917 else { prop = 0.; }
2918
2919 // Time walk correction
2920 double Tw = 0;
2921 // if(expmc==1) {
2922 // Calmdc_const3 &TwRef = Cal3Mgr[layer];
2923 // if(adc[wid]>0.) Tw = TwRef.tw(0) + TwRef.tw(1)/sqrt(adc[wid]);
2924 // }
2925
2926 double driftt = 0;
2927 double dummy;
2928
2929 double xw = fx + ( bx - fx ) / ( bz - fz ) * ( z - fz );
2930 double yw = fy + ( by - fy ) / ( bz - fz ) * ( z - fz );
2931 // move pivot to the wire (slant angle ignored)
2932 HepPoint3D pivot1( xw, yw, z );
2933 helix0.pivot( pivot1 );
2934
2935 double zw = helix0.a()[3]; // z position
2936
2937 int lr = 2;
2938 // if((xw*helix0.x(0.).y()-yw*helix0.x(0.).x())<0) lr=-1;
2939 double p[3];
2940 p[0] = helix0.momentum( 0. ).x();
2941 p[1] = helix0.momentum( 0. ).y();
2942 double pos[2];
2943 pos[0] = xw;
2944 pos[1] = yw;
2945 double alpha;
2946 //// calmdc_getalpha_( p, pos, &alpha);
2947 // double dist = wir.distance(); this is wrong
2948 // double dist = abs(helix0.a()[0]); this is wrong
2949 double dist = 0.;
2950
2951 dist = ( m_mdcUtilitySvc->doca( layer, wid, helix_par, helixErr ) ) * 10.0;
2952
2953 if ( dist < 0. ) lr = 1;
2954 else lr = 0;
2955 dist = fabs( dist );
2956 if ( dist > 0.4 * ( mdcGeomSvc->Layer( layer ) )->PCSiz() ) continue;
2957
2958 //* drift distance cut
2959 // int ia;
2960 // ia = (int) ((alpha+90.)/10.);
2961 // double da = alpha+90. - ia*10.;
2962 // if (ia==0) ia=18;
2963
2964 // Calmdc_const &BndRef = Cal1Mgr[layer];
2965 // if(lr==-1) {
2966 // if(dist < BndRef.bndl(ia,0)) continue;
2967 // if(dist > BndRef.bndl(ia,3)) continue;
2968 // }
2969 // if(lr== 1) {
2970 // if(dist < BndRef.bndr(ia,0)) continue;
2971 // if(dist > BndRef.bndr(ia,3)) continue;
2972 // }
2973
2974 if ( !m_useXT ) { driftt = dist / ( Estparam.vdrift() ); }
2975 else
2976 {
2977 double entrance = ( *it_gothit )->getEntra();
2978 driftt = m_mdcCalibFunSvc->distToDriftTime( dist, layer, wid, lr, entrance );
2979 }
2980 if ( m_useT0cal )
2981 {
2982 T0_cal = m_mdcCalibFunSvc->getT0( layer, wid ) +
2983 m_mdcCalibFunSvc->getTimeWalk( layer, tdc );
2984 }
2985
2986 double zprop = fabs( zw - m_zst[layer] );
2987 double tp = zprop / m_vp[layer];
2988 // cout<<"proptation time --tp st: "<<tp<<endl;
2989 if ( driftt > tdc ) continue;
2990 double difft = tdc - driftt - toft - tp - T0_cal;
2991 if ( difft < -10 ) continue;
2992 if ( ndriftt >= 500 ) break;
2993 Mdc_t0Arr[ndriftt] = difft;
2994 // if(difft>-2 && difft<22)
2995 // if(fabs(hitChi2)<0.2)
2996 // if(difft>-20 && difft<30)
2997 sum_EstimeMdc = sum_EstimeMdc + difft;
2998 ndriftt++;
2999
3000 // htf[2]->accumulate( t0_minus_TDC[wid] );
3001
3002 double tev = -t0_minus_TDC[wid] + driftt;
3003 if ( Estparam.MDC_Tof() != 0 ) tev += direction * toft;
3004 if ( Estparam.MDC_Prop() != 0 ) tev += prop;
3005 //// if(Estparam.MDC_Walk()!=0) tev += Tw;
3006
3007 // sum_tev_st += tev;
3008
3009 // htf[3+hid]->accumulate( tev );
3010
3011 nhits_st++;
3012 tev_st[nhits_st - 1] = tev;
3013
3014 if ( Estparam.MDC_Debug() != 0 )
3015 log << MSG::INFO << "*** tev_st ***" << tev << endmsg;
3016 double driftt_mea = t0_minus_TDC[wid];
3017 // if(abs(driftt-t0_minus_TDC[wid]) <75.)
3018 if ( abs( driftt - driftt_mea ) < 75. )
3019 {
3020 // sum_tev_st2 += tev;
3021 nhits_st2++;
3022 if ( Estparam.MDC_Debug() != 0 )
3023 log << MSG::INFO << "*** tev_st2 ***" << tev << endmsg;
3024 }
3025 } //* End stereo hits
3026
3027 } // End hits
3028 nmatch_mdc++;
3029 } //* End tracks
3030
3031 if ( m_ntupleflag && m_tuple2 ) g_nmatchmdc = nmatch_mdc;
3032 if ( ndriftt != 0 )
3033 {
3034 if ( m_mdcopt ) { sum_EstimeMdc = Opt_new( Mdc_t0Arr, ndriftt, 400.0 ); }
3035 else { sum_EstimeMdc = sum_EstimeMdc / ndriftt; }
3036 if ( m_ntupleflag && m_tuple2 ) g_EstimeMdc = sum_EstimeMdc;
3037 t_Est = sum_EstimeMdc + tOffset_b; // 12.yzhang add tOffset_b, calc. Tes after rec for ?
3038 if ( t_Est < 0 ) t_Est = 0;
3039 if ( optCosmic == 0 )
3040 {
3041 tEstFlag = 102;
3042 nbunch = ( (int)( t_Est - offset ) ) / bunchtime;
3043 // if(((int)(t_Est-offset))%bunchtime>bunchtime/2) nbunch=nbunch+1;
3044 if ( ( t_Est - offset - nbunch * bunchtime ) > ( bunchtime / 2 ) ) nbunch = nbunch + 1;
3045 t_Est = nbunch * bunchtime + offset + tOffset_b;
3046 // tEstFlag=2;
3047 /* if(m_nbunch==6){
3048 nbunch=((int)(t_Est-offset))/4;
3049 if(((int)(t_Est-offset))%4>2 ||((int)(t_Est-offset))%4==2 ) nbunch=nbunch+1;
3050 t_Est=nbunch*4+offset;
3051 tEstFlag=2;
3052 }*/
3053 }
3054 if ( optCosmic )
3055 {
3056 t_Est = sum_EstimeMdc;
3057 tEstFlag = 202;
3058 }
3059 }
3060 if ( m_ntupleflag && m_tuple2 ) g_ndriftt = ndriftt;
3061 }
3062 // yzhang add, Store to TDS
3063 if ( t_Est != -999 )
3064 {
3065 // changed event start time flag after rec
3066 if ( ( !m_beforrec ) && ( Testime_fzisan != t_Est ) )
3067 {
3068 if ( tEstFlag == 211 ) tEstFlag = 213;
3069 if ( tEstFlag == 212 ) tEstFlag = 216;
3070 if ( tEstFlag == 111 ) tEstFlag = 113;
3071 if ( tEstFlag == 112 ) tEstFlag = 116;
3072 }
3073
3074 if ( optCosmic /*&& (nmatch_barrel || nmatch_end)*/ )
3075 {
3076 StatusCode scStoreTds = storeTDS( t_Est, tEstFlag, t_quality );
3077 if ( scStoreTds != StatusCode::SUCCESS ) { return scStoreTds; }
3078 }
3079 else if ( !optCosmic )
3080 {
3081 StatusCode scStoreTds = storeTDS( t_Est, tEstFlag, t_quality );
3082 if ( scStoreTds != StatusCode::SUCCESS ) { return scStoreTds; }
3083 }
3084 }
3085 else
3086 {
3087 // no t_Est calc.
3088 if ( m_beforrec )
3089 {
3090 // store event start time from segment fitting method
3091 // FIXME add tOffset_b or tOffset_e ???
3092 double segTest = Testime_fzisan + tOffset_b;
3093 int segFlag = TestimeFlag_fzisan;
3094 double segQuality = TestimeQuality_fzisan;
3095 StatusCode scStoreTds = storeTDS( segTest, segFlag, segQuality );
3096 if ( scStoreTds != StatusCode::SUCCESS ) { return scStoreTds; }
3097 }
3098 }
3099 // zhangy add, end of Store to TDS
3100
3101 // check RecEsTimeCol registered
3102 SmartDataPtr<RecEsTimeCol> arecestimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
3103 if ( !arecestimeCol )
3104 {
3105 if ( m_debug == 4 ) log << MSG::WARNING << "Could not find RecEsTimeCol" << endmsg;
3106 return ( StatusCode::SUCCESS );
3107 }
3108 RecEsTimeCol::iterator iter_evt = arecestimeCol->begin();
3109 for ( ; iter_evt != arecestimeCol->end(); iter_evt++ )
3110 {
3111 log << MSG::INFO << "Test = " << ( *iter_evt )->getTest()
3112 << ", Status = " << ( *iter_evt )->getStat() << endmsg;
3113 if ( m_ntupleflag && m_tuple2 ) { g_Testime = ( *iter_evt )->getTest(); }
3114 // cout<<"register to TDS: "<<(*iter_evt)->getTest()<<", "<<(*iter_evt)->getStat()<<endl;
3115 }
3116
3117 if ( m_ntupleflag )
3118 {
3119 if ( m_tuple2 )
3120 {
3121 for ( g_ntofdown = 0; g_ntofdown < ntofdown; g_ntofdown++ )
3122 { g_meantevdown[g_ntofdown] = meantevdown[g_ntofdown]; }
3123 for ( g_ntofup = 0; g_ntofup < ntofup; g_ntofup++ )
3124 { g_meantevup[g_ntofup] = meantevup[g_ntofup]; }
3125 g_nmatch_tot = nmatch;
3126 m_estFlag = tEstFlag; /////20100427 guanyh add
3127 m_estTime = t_Est;
3128 StatusCode status = m_tuple2->write();
3129 if ( !status.isSuccess() ) { log << MSG::ERROR << "Can't fill ntuple!" << endmsg; }
3130 }
3131 if ( m_tuple9 )
3132 {
3133 for ( g_nmatch = 0; g_nmatch < nmatch; g_nmatch++ )
3134 { g_meantev[g_nmatch] = meantev[g_nmatch]; }
3135 StatusCode status = m_tuple9->write();
3136 if ( !status.isSuccess() ) { log << MSG::ERROR << "Can't fill ntuple!" << endmsg; }
3137 }
3138 }
3139 return StatusCode::SUCCESS;
3140}
3141
3142StatusCode EsTimeAlg::finalize() {
3143
3144 MsgStream log( msgSvc(), name() );
3145 log << MSG::INFO << "in finalize()" << endmsg;
3146 if ( m_ntupleflag && m_tuple3 )
3147 {
3148 StatusCode status = m_tuple3->write();
3149 if ( !status.isSuccess() ) { log << MSG::ERROR << "Can't fill ntuple!" << endmsg; }
3150 }
3151 cout << "EsTimeAlg::finalize(),total events in this run: " << m_pass[0] << endl;
3152 return StatusCode::SUCCESS;
3153}
3154
3155// make TDS
3156StatusCode EsTimeAlg::storeTDS( double tEst, int tEstFlag, double quality ) {
3157 StatusCode sc;
3158 MsgStream log( msgSvc(), name() );
3159
3160 // check whether Recon already registered
3161 DataObject* aReconEvent;
3162 eventSvc()->findObject( "/Event/Recon", aReconEvent );
3163 if ( aReconEvent == NULL )
3164 {
3165 // then register Recon
3166 aReconEvent = new ReconEvent();
3167 sc = eventSvc()->registerObject( "/Event/Recon", aReconEvent );
3168 if ( sc != StatusCode::SUCCESS )
3169 {
3170 log << MSG::FATAL << "Could not register ReconEvent" << endmsg;
3171 return StatusCode::FAILURE;
3172 }
3173 }
3174
3175 // get the DataManager service
3176 SmartIF<IDataManagerSvc> dataManagerSvc( eventSvc() );
3177
3178 // clear MdcFastTrk
3179 DataObject* aRecMdcTrack;
3180 eventSvc()->findObject( "/Event/Recon/RecMdcTrackCol", aRecMdcTrack );
3181 if ( aRecMdcTrack != NULL )
3182 {
3183 dataManagerSvc->clearSubTree( "/Event/Recon/RecMdcTrackCol" );
3184 aRecMdcTrack = NULL;
3185 }
3186
3187 if ( tEst < 0 ) { return StatusCode::SUCCESS; }
3188
3189 // clear RecEsTimeCol
3190 DataObject* aRecEsTime;
3191 eventSvc()->findObject( "/Event/Recon/RecEsTimeCol", aRecEsTime );
3192 if ( aRecEsTime != NULL )
3193 {
3194 dataManagerSvc->clearSubTree( "/Event/Recon/RecEsTimeCol" );
3195 aRecEsTime = NULL;
3196 }
3197
3198 // register event start time to TDS
3199 RecEsTimeCol* aRecEsTimeCol = new RecEsTimeCol;
3200 sc = eventSvc()->registerObject( "/Event/Recon/RecEsTimeCol", aRecEsTimeCol );
3201 if ( sc != StatusCode::SUCCESS )
3202 {
3203 log << MSG::ERROR << "Could not register RecEsTimeCol" << endmsg;
3204 return StatusCode::FAILURE;
3205 }
3206
3207 RecEsTime* arecestime = new RecEsTime;
3208 arecestime->setTest( tEst );
3209 arecestime->setStat( tEstFlag );
3210 arecestime->setQuality( quality );
3211 aRecEsTimeCol->push_back( arecestime );
3212
3213 return StatusCode::SUCCESS;
3214}
3215
3216double EsTimeAlg::Opt_new( const double* arr, const int size, const double sigma_cut ) {
3217 Vdouble t0v_mdc;
3218 t0v_mdc.clear();
3219 double mean = -999;
3220 double sigma = 9999;
3221 for ( int i = 0; i < size; i++ ) { t0v_mdc.push_back( arr[i] ); }
3222 if ( size == 0 ) mean = -999;
3223 if ( size == 1 ) mean = t0v_mdc[0];
3224 if ( size == 2 ) mean = 0.5 * ( t0v_mdc[0] + t0v_mdc[1] );
3225 if ( size >= 3 )
3226 {
3227 for ( int n = 0; n < size; n++ )
3228 {
3229 mean = 0.0;
3230 sigma = 0.0;
3231 for ( int i = 0; i < t0v_mdc.size(); i++ ) { mean += t0v_mdc[i]; }
3232 mean = mean / t0v_mdc.size();
3233 for ( int i = 0; i < t0v_mdc.size(); i++ )
3234 { sigma += ( t0v_mdc[i] - mean ) * ( t0v_mdc[i] - mean ); }
3235 sigma = sigma / t0v_mdc.size();
3236 if ( sigma < sigma_cut ) break;
3237 double tmp = fabs( mean - t0v_mdc[0] );
3238 int no = 0;
3239 for ( int j = 0; j < t0v_mdc.size(); j++ )
3240 {
3241 if ( fabs( mean - t0v_mdc[j] ) >= tmp )
3242 {
3243 no = j;
3244 tmp = fabs( mean - t0v_mdc[j] );
3245 }
3246 }
3247 t0v_mdc.erase( t0v_mdc.begin() + no );
3248 if ( t0v_mdc.size() <= 2 ) break;
3249 }
3250 mean = 0.0;
3251 for ( int i = 0; i < t0v_mdc.size(); i++ ) { mean += t0v_mdc[i]; }
3252 mean = mean / t0v_mdc.size();
3253 }
3254 return mean;
3255}
3256
3257double EsTimeAlg::EST_Trimmer( double t0_original, double t0_offset, double raw_offset,
3258 double t0_offset_dt, double bunchTime ) {
3259 int Nbunch = (int)( t0_original - t0_offset - raw_offset ) / bunchTime;
3260 if ( ( t0_original - t0_offset - raw_offset - bunchTime * Nbunch ) > ( bunchTime / 2. ) )
3261 { Nbunch = Nbunch + 1; }
3262 double t_Est = Nbunch * bunchTime + t0_offset + t0_offset_dt;
3263 return t_Est;
3264}
momentum
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
Definition BesBdkRc/src/fortran/HepEvt.h:23
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
runNo
int runNo
Definition DQA_TO_DB.cxx:13
n
const Int_t n
Definition DataBase/tau_mode.c:57
difft
Double_t difft
Definition DataBase/tau_mode.c:34
Emc_helix.h
VELPROP
const double VELPROP
Definition EsTimeAlg.cxx:81
PROTONMAS2
const double PROTONMAS2
Definition EsTimeAlg.cxx:76
MXTRK
#define MXTRK
MUMAS2
const double MUMAS2
Definition EsTimeAlg.cxx:74
PIMAS2
const double PIMAS2
Definition EsTimeAlg.cxx:75
ELMAS2
const double ELMAS2
Definition EsTimeAlg.cxx:73
VLIGHT
const double VLIGHT
Definition EsTimeAlg.cxx:72
RCTOF2
const double RCTOF2
Definition EsTimeAlg.cxx:77
TDC2NSEC
const double TDC2NSEC
Definition EsTimeAlg.cxx:79
MXTKHIT
#define MXTKHIT
MXWIRE
#define MXWIRE
PIBY44
const double PIBY44
Definition EsTimeAlg.cxx:80
RCEMC2
const double RCEMC2
Definition EsTimeAlg.cxx:78
NSEC2TDC
const double NSEC2TDC
Definition EsTimeAlg.cxx:79
Vdouble
std::vector< double > Vdouble
Definition EsTimeAlg.h:25
RecEsTimeCol
ObjectVector< RecEsTime > RecEsTimeCol
Definition Event/EvTimeEvent/include/EvTimeEvent/RecEsTime.h:48
PI
#define PI
Definition EvtD0ToKKpi0.cc:40
alpha
double alpha
Definition EvtPhokharaDef.hh:17
Vdouble
std::vector< double > Vdouble
Definition Gam4pikp.cxx:39
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcRecEvent/RecMdcTrack.h:22
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition KarFin.h:34
eventNo
int eventNo
Definition MdcDedxRecon.cxx:50
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
M_PI
#define M_PI
Definition TConstant.h:4
detVerSvc
IDetVerSvc * detVerSvc
Definition TofRawDataProvider.cxx:30
Toffz_helix.h
Emc_helix
Definition Emc_helix.h:15
Emc_helix::Z_emc
double Z_emc
Definition Emc_helix.h:28
Emc_helix::theta_emc
double theta_emc
Definition Emc_helix.h:29
Emc_helix::phi_emc
double phi_emc
Definition Emc_helix.h:30
Emc_helix::Emc_Get
int Emc_Get(double, int, double[])
Definition Emc_helix.cxx:65
Emc_helix::pathlCut
void pathlCut(double pathl_max)
Definition Emc_helix.h:42
EsTimeAlg::offset_dt
double offset_dt
Definition EsTimeAlg.h:54
EsTimeAlg::m_useXT
bool m_useXT
Definition EsTimeAlg.h:78
EsTimeAlg::m_optCosmic
int m_optCosmic
Definition EsTimeAlg.h:46
EsTimeAlg::m_ntupleflag
int m_ntupleflag
Definition EsTimeAlg.h:43
EsTimeAlg::m_useTimeFactor
bool m_useTimeFactor
Definition EsTimeAlg.h:84
EsTimeAlg::m_cosmicScheme
int m_cosmicScheme
Definition EsTimeAlg.h:47
EsTimeAlg::m_useSw
bool m_useSw
Definition EsTimeAlg.h:80
EsTimeAlg::m_bunchtime_MC
double m_bunchtime_MC
Definition EsTimeAlg.h:42
EsTimeAlg::m_TofOpt_Cut
double m_TofOpt_Cut
Definition EsTimeAlg.h:83
EsTimeAlg::m_beforrec
int m_beforrec
Definition EsTimeAlg.h:57
EsTimeAlg::toffset_rawtime_e
double toffset_rawtime_e
Definition EsTimeAlg.h:53
EsTimeAlg::EsTimeAlg
EsTimeAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition EsTimeAlg.cxx:83
EsTimeAlg::m_ebeam
double m_ebeam
Definition EsTimeAlg.h:61
EsTimeAlg::offset_dt_end
double offset_dt_end
Definition EsTimeAlg.h:55
EsTimeAlg::m_userawtime_opt
int m_userawtime_opt
Definition EsTimeAlg.h:51
EsTimeAlg::m_mdcopt
bool m_mdcopt
Definition EsTimeAlg.h:81
EsTimeAlg::finalize
StatusCode finalize()
Definition EsTimeAlg.cxx:3142
EsTimeAlg::m_TofOpt
bool m_TofOpt
Definition EsTimeAlg.h:82
EsTimeAlg::m_phase
int m_phase
Definition EsTimeAlg.h:50
EsTimeAlg::initialize
StatusCode initialize()
Definition EsTimeAlg.cxx:122
EsTimeAlg::m_evtNo
int m_evtNo
Definition EsTimeAlg.h:60
EsTimeAlg::m_flag
int m_flag
Definition EsTimeAlg.h:37
EsTimeAlg::m_nbunch
int m_nbunch
Definition EsTimeAlg.h:40
EsTimeAlg::toffset_rawtime
double toffset_rawtime
Definition EsTimeAlg.h:52
EsTimeAlg::m_useT0cal
bool m_useT0cal
Definition EsTimeAlg.h:79
EsTimeAlg::execute
StatusCode execute()
Definition EsTimeAlg.cxx:308
EsTimeAlg::m_debug
int m_debug
Definition EsTimeAlg.h:56
EstParameter::pathlCut
double pathlCut() const
Definition EstParameter.cxx:83
EstParameter::ztofCutmin
double ztofCutmin() const
Definition EstParameter.cxx:85
EstParameter::MDC_Prop
double MDC_Prop() const
Definition EstParameter.cxx:101
EstParameter::MDC_drCut
double MDC_drCut() const
Definition EstParameter.cxx:95
EstParameter::ztofCutmax
double ztofCutmax() const
Definition EstParameter.cxx:87
EstParameter::MDC_dzCut
double MDC_dzCut() const
Definition EstParameter.cxx:97
EstParameter::MDC_Debug
int MDC_Debug() const
Definition EstParameter.cxx:105
EstParameter::MDC_Tof
double MDC_Tof() const
Definition EstParameter.cxx:99
EstParameter::vdrift
double vdrift() const
Definition EstParameter.cxx:71
EstParameter::MDC_Inner
double MDC_Inner() const
Definition EstParameter.cxx:93
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:197
Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition Reconstruction/TrackUtil/src/Helix.cxx:172
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Reconstruction/TrackUtil/src/Helix.cxx:120
Helix::radius
double radius(void) const
returns radious of helix.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:201
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:217
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:199
IMdcGeomSvc::Layer
virtual const MdcGeoLayer *const Layer(unsigned id)=0
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
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
RawDataUtil::MdcTime
static double MdcTime(int timeChannel)
Definition Event/RawEvent/include/RawEvent/RawDataUtil.h:7
RawDataUtil::TofTime
static double TofTime(unsigned int timeChannel)
Definition Event/RawEvent/include/RawEvent/RawDataUtil.h:57
TofFz_helix
Definition Toffz_helix.h:17
TofFz_helix::TofFz_Get
int TofFz_Get(double, int, double[])
Definition Toffz_helix.cxx:70
TofFz_helix::Kappa
double Kappa
Definition Toffz_helix.h:38
TofFz_helix::r_endtof
double r_endtof
Definition Toffz_helix.h:41
TofFz_helix::Pathl
double Pathl
Definition Toffz_helix.h:31
TofFz_helix::Path_etf
double Path_etf
Definition Toffz_helix.h:33
TofFz_helix::r_etf
double r_etf
Definition Toffz_helix.h:42
TofFz_helix::Tofid
int Tofid
Definition Toffz_helix.h:26
TofFz_helix::Etfid
int Etfid
Definition Toffz_helix.h:27
TofFz_helix::ztofCut
void ztofCut(double ztof_min, double ztof_max)
Definition Toffz_helix.h:51
TofFz_helix::pathlCut
void pathlCut(double pathl_max)
Definition Toffz_helix.h:49
TofFz_helix::Tanl
double Tanl
Definition Toffz_helix.h:38
TofFz_helix::Z_tof
double Z_tof
Definition Toffz_helix.h:34
TofFz_helix::Z_etf
double Z_etf
Definition Toffz_helix.h:35
TofID::endcap
static int endcap(const Identifier &id)
Definition TofID.cxx:108
t
int t()
Definition t.c:1