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