CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemMdcCombAlg.cxx
Go to the documentation of this file.
1#include "CgemMdcCombAlg/CgemMdcCombAlg.h"
2#include "GaudiKernel/IDataManagerSvc.h"
3#include "GaudiKernel/SmartDataPtr.h"
4#include "GaudiKernel/IDataProviderSvc.h"
5#include "GaudiKernel/IPartPropSvc.h"
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/AlgFactory.h"
8#include "GaudiKernel/ISvcLocator.h"
9#include "GaudiKernel/DataSvc.h"
10#include "GaudiKernel/PropertyMgr.h"
11#include "GaudiKernel/Bootstrap.h"
12#include "GaudiKernel/INTupleSvc.h"
13
14
15#include "EventModel/EventHeader.h"
16#include "CLHEP/Units/PhysicalConstants.h"
17
18#include "HepPDT/ParticleDataTable.hh"
19#include "HepPDT/ParticleData.hh"
20
21#include "EventModel/Event.h"
22#include "EvTimeEvent/RecEsTime.h"
23#include "EvtRecEvent/EvtRecEvent.h"
24#include "EvtRecEvent/EvtRecTrack.h"
25#include "Identifier/Identifier.h"
26#include "RawDataProviderSvc/IRawDataProviderSvc.h"
27#include "RawDataProviderSvc/MdcRawDataProvider.h"
28#include "RawDataProviderSvc/RawDataProviderSvc.h"
29#include "CLHEP/Vector/ThreeVector.h"
30#include "CLHEP/Matrix/Vector.h"
31#include "CLHEP/Matrix/SymMatrix.h"
32#include "CLHEP/Geometry/Point3D.h"
33#include "CLHEP/Vector/LorentzVector.h"
34//#include "TrackUtil/Helix.h"
35//#include "TrkFitter/TrkHelixFitter.h"
36//#include "TrkFitter/TrkHelixMaker.h"
37#include "KalFitAlg/helix/Helix.h"
38
39#include "CLHEP/Matrix/Matrix.h"
40
41#include "CgemGeomSvc/CgemGeomSvc.h"
42#include "CgemRecEvent/RecCgemCluster.h"
43#include "CgemRecEvent/RecCgemSegment.h"
44#include "CgemMdcCombAlg/CgemMdcCombAlg.h"
45#include "McTruth/CgemMcHit.h"
46
47#include "MdcRecEvent/RecMdcTrack.h"
48#include "MdcRecEvent/RecMdcHit.h"
49
50#include "TRandom.h"
51#include "TArrayI.h"
52#include "TF1.h"
53#include "TMath.h"
54#include "TMinuit.h"
55
56
57#include <vector>
58#include <iostream>
59#include <fstream>
60#include <iomanip>
61#include <cstring>
62
63using namespace std;
64using namespace Event;
65
66//////////////////////////////////////////////////////////////////////////////////
67CgemMdcCombAlg::CgemMdcCombAlg(const std::string& name, ISvcLocator* pSvcLocator) :
68 Algorithm(name, pSvcLocator)
69{
70 declareProperty("debug", m_debug = 0);
71 declareProperty("prodNT", m_prodNT = 0);
72 declareProperty("dropHot", m_dropHot= 1);
73 declareProperty("combineTracking",m_combineTracking=0);
74 declareProperty("keepBadTdc",m_keepBadTdc=0);
75 declareProperty("keepUnmatch",m_keepUnmatch=1);
76 declareProperty("dr_par",m_dr_par);
77 declareProperty("phi0_par",m_phi0_par);
78 declareProperty("kappa_par",m_kappa_par);
79 declareProperty("dz_par",m_dz_par);
80 declareProperty("tanl_par",m_tanl_par);
81 declareProperty("Pt",m_Pt);
82 declareProperty("w_dr",m_w_dr);
83 declareProperty("w_phi0",m_w_phi0);
84 declareProperty("w_kappa",m_w_kappa);
85 declareProperty("w_dz",m_w_dz);
86 declareProperty("w_tanl",m_w_tanl);
87
88 declareProperty("P_T", mv_Pt);
89 declareProperty("mean_dr", mv_mean_dr);
90 declareProperty("mean_phi0", mv_mean_phi0);
91 declareProperty("mean_kappa",mv_mean_kappa);
92 declareProperty("mean_dz", mv_mean_dz);
93 declareProperty("mean_tanL", mv_mean_tanL);
94 declareProperty("sig_dr", mv_sig_dr);
95 declareProperty("sig_phi0", mv_sig_phi0);
96 declareProperty("sig_kappa", mv_sig_kappa);
97 declareProperty("sig_dz", mv_sig_dz);
98 declareProperty("sig_tanL", mv_sig_tanL);
99
100 declareProperty("matchChi2Cut", mv_matchChi2Cut);
101}
102
103// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
104StatusCode CgemMdcCombAlg::initialize()
105{
106
107 MsgStream log(msgSvc(), name());
108 log << MSG::INFO << "in initialize()" << endreq;
109 if(m_prodNT)
110 {
111 /////////////////////////////// book ntuple //////////////////////////////////////
112 INTupleSvc* ntupleSvc;
113 StatusCode sc = service("NTupleSvc", ntupleSvc);
114 NTuplePtr mt(ntupleSvc, "FILE199/CgemMdcComb");
115 if ( mt ) match_ntuple = mt;
116 else{
117 match_ntuple = ntupleSvc->book("FILE199/CgemMdcComb", CLID_ColumnWiseTuple, "match");
118 StatusCode status;
119 if(match_ntuple){
120 status = match_ntuple->addItem("run", run);
121 status = match_ntuple->addItem("evt", evt);
122 status = match_ntuple->addItem("costhetaMC", costhetaMC);
123 status = match_ntuple->addItem("phiMC", phiMC);
124 status = match_ntuple->addItem("pdgMC", pdgMC);
125 status = match_ntuple->addItem("truth_dr", truth_dr);
126 status = match_ntuple->addItem("truth_dz", truth_dz);
127 status = match_ntuple->addItem("truth_phi0", truth_phi0);
128 status = match_ntuple->addItem("truth_kappa", truth_kappa);
129 status = match_ntuple->addItem("truth_tanl", truth_tanl);
130 status = match_ntuple->addItem("track_size",track_size,0,1000);
131 status = match_ntuple->addIndexedItem("matched", track_size, CgemSegMatched);
132 status = match_ntuple->addIndexedItem("dr_avg", track_size, dr_avg);
133 status = match_ntuple->addIndexedItem("dz_avg", track_size, dz_avg);
134 status = match_ntuple->addIndexedItem("tanL_avg", track_size, tanL_avg);
135 status = match_ntuple->addItem("segment_size",segment_size,0,1000);
136 status = match_ntuple->addIndexedItem("fit_failed", segment_size , match_failed);
137 status = match_ntuple->addItem("nMatch",match_nsegment,0,1000);
138 status = match_ntuple->addIndexedItem("chisq_dr", match_nsegment , match_dr_chisq);
139 status = match_ntuple->addIndexedItem("chisq_phi0", match_nsegment , match_phi0_chisq);
140 status = match_ntuple->addIndexedItem("chisq_kappa", match_nsegment , match_kappa_chisq);
141 status = match_ntuple->addIndexedItem("chisq_tanl", match_nsegment , match_tanl_chisq);
142 status = match_ntuple->addIndexedItem("chisq_dz", match_nsegment , match_dz_chisq);
143 status = match_ntuple->addIndexedItem("IsTruth", match_nsegment , match_IsTruth);
144 status = match_ntuple->addIndexedItem("track_dr", match_nsegment , track_dr);
145 status = match_ntuple->addIndexedItem("track_phi0", match_nsegment , track_phi0);
146 status = match_ntuple->addIndexedItem("track_kappa", match_nsegment , track_kappa);
147 status = match_ntuple->addIndexedItem("track_tanl", match_nsegment , track_tanl);
148 status = match_ntuple->addIndexedItem("track_dz", match_nsegment , track_dz);
149 status = match_ntuple->addIndexedItem("track_dr_err", match_nsegment , track_dr_err);
150 status = match_ntuple->addIndexedItem("track_phi0_err", match_nsegment , track_phi0_err);
151 status = match_ntuple->addIndexedItem("track_kappa_err", match_nsegment , track_kappa_err);
152 status = match_ntuple->addIndexedItem("track_tanl_err", match_nsegment , track_tanl_err);
153 status = match_ntuple->addIndexedItem("track_dz_err", match_nsegment , track_dz_err);
154 status = match_ntuple->addIndexedItem("seg_dr", match_nsegment , seg_dr);
155 status = match_ntuple->addIndexedItem("seg_phi0", match_nsegment , seg_phi0);
156 status = match_ntuple->addIndexedItem("seg_kappa", match_nsegment , seg_kappa);
157 status = match_ntuple->addIndexedItem("seg_tanl", match_nsegment , seg_tanl);
158 status = match_ntuple->addIndexedItem("seg_dz", match_nsegment , seg_dz);
159 status = match_ntuple->addIndexedItem("seg_dr_err", match_nsegment , seg_dr_err);
160 status = match_ntuple->addIndexedItem("seg_phi0_err", match_nsegment , seg_phi0_err);
161 status = match_ntuple->addIndexedItem("seg_kappa_err", match_nsegment , seg_kappa_err);
162 status = match_ntuple->addIndexedItem("seg_tanl_err", match_nsegment , seg_tanl_err);
163 status = match_ntuple->addIndexedItem("seg_dz_err", match_nsegment , seg_dz_err);
164 }
165 }
166 }
167 unSegNum = 0;
168 unTrkNum = 0;
169 EvInWindow = 0;
170 EvOutWindow = 0;
171
172 mg_mean_dr=0;
173 mg_mean_phi0=0;
174 mg_mean_kappa=0;
175 mg_mean_dz=0;
176 mg_mean_tanL=0;
177 mg_sig_dr=0;
178 mg_sig_phi0=0;
179 mg_sig_kappa=0;
180 mg_sig_dz=0;
181 mg_sig_tanL=0;
182
183 m_nPt=0;
184 if(mv_Pt.size()>0) {
185 cout<<"mv_Pt.size="<<mv_Pt.size()<<endl;
186 m_nPt = mv_Pt.size();
187 cout<<"pt[0]="<<mv_Pt[0]<<", pt[npt-1]="<<mv_Pt[m_nPt-1]<<endl;
188 if(mv_mean_dr.size()==m_nPt) mg_mean_dr = new TGraph(m_nPt, &mv_Pt[0], &mv_mean_dr[0]);
189 else cout<<"size of mean_dr != size of PT !"<<endl;
190
191 if(mv_mean_phi0.size()==m_nPt) mg_mean_phi0 = new TGraph(m_nPt, &mv_Pt[0], &mv_mean_phi0[0]);
192 else cout<<"size of mean_phi0 != size of PT !"<<endl;
193
194 if(mv_mean_kappa.size()==m_nPt)mg_mean_kappa= new TGraph(m_nPt, &mv_Pt[0], &mv_mean_kappa[0]);
195 else cout<<"size of mean_kappa != size of PT !"<<endl;
196
197 if(mv_mean_dz.size()==m_nPt) mg_mean_dz = new TGraph(m_nPt, &mv_Pt[0], &mv_mean_dz[0]);
198 else cout<<"size of mean_dz != size of PT !"<<endl;
199
200 if(mv_mean_tanL.size()==m_nPt) mg_mean_tanL = new TGraph(m_nPt, &mv_Pt[0], &mv_mean_tanL[0]);
201 else cout<<"size of mean_tanL != size of PT !"<<endl;
202
203 if(mv_sig_dr.size()==m_nPt) mg_sig_dr = new TGraph(m_nPt, &mv_Pt[0], &mv_sig_dr[0]);
204 else cout<<"size of sig_dr != size of PT !"<<endl;
205
206 if(mv_sig_phi0.size()==m_nPt) mg_sig_phi0 = new TGraph(m_nPt, &mv_Pt[0], &mv_sig_phi0[0]);
207 else cout<<"size of sig_phi0 != size of PT !"<<endl;
208
209 if(mv_sig_kappa.size()==m_nPt) mg_sig_kappa = new TGraph(m_nPt, &mv_Pt[0], &mv_sig_kappa[0]);
210 else cout<<"size of sig_kappa != size of PT !"<<endl;
211
212 if(mv_sig_dz.size()==m_nPt) mg_sig_dz = new TGraph(m_nPt, &mv_Pt[0], &mv_sig_dz[0]);
213 else cout<<"size of sig_dz != size of PT !"<<endl;
214
215 if(mv_sig_tanL.size()==m_nPt) mg_sig_tanL = new TGraph(m_nPt, &mv_Pt[0], &mv_sig_tanL[0]);
216 else cout<<"size of sig_tanL != size of PT !"<<endl;
217 }
218
219 return StatusCode::SUCCESS;
220
221}
222
223StatusCode CgemMdcCombAlg::execute()
224{
225 exe_v2();
226
227 return StatusCode::SUCCESS;
228}
229
230void CgemMdcCombAlg::getMcPar() {
231 MsgStream log(msgSvc(), name());
232 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
233 if (!mcParticleCol) {
234 log << MSG::ERROR<< "Could not find McParticle" << endreq;
235 } else{
236 if(m_debug==1) cout<< "size_mcParticleCol "<<mcParticleCol->size()<< endl;
237 McParticleCol::iterator iter_mc = mcParticleCol->begin();
238 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
239 bool findJpsi=false;
240 int pdg=(*iter_mc)->particleProperty();
241 if(abs(pdg)==13){
242 findJpsi=true;
243 }
244 if(!findJpsi) continue;
245 double mc_charge;
246
247 if(pdg<0)mc_charge = 1.;
248 else mc_charge =-1.;
249 double truth_x = (*iter_mc)->initialPosition().x();
250 double truth_y = (*iter_mc)->initialPosition().y();
251 double truth_z = (*iter_mc)->initialPosition().z();
252
253 double px = (*iter_mc)->initialFourMomentum().px();
254 double py = (*iter_mc)->initialFourMomentum().py();
255 double pz = (*iter_mc)->initialFourMomentum().pz();
256
257 HepPoint3D truth_position(truth_x, truth_y, truth_z);
258 Hep3Vector p3(px, py, pz);
259 KalmanFit::Helix* tmp_helix = new KalmanFit::Helix(truth_position, p3, mc_charge);
260 HepPoint3D tmp_pivot(0., 0., 0.);
261 tmp_helix->pivot(tmp_pivot);
262
263 if(m_prodNT) {
264 truth_dr = tmp_helix->dr();
265 truth_dz = tmp_helix->dz();
266 truth_phi0 = tmp_helix->phi0();
267 truth_kappa = tmp_helix->kappa();
268 truth_tanl = tmp_helix->tanl();
269
270 costhetaMC = pz/sqrt(pow(px,2)+pow(py,2)+pow(pz,2));
271 phiMC = p3.phi();
272 pdgMC = pdg;
273 }
274 delete tmp_helix;
275 }
276 }
277}//getMcPar()
278
279StatusCode CgemMdcCombAlg::exe_v1()
280{
281 MsgStream log(msgSvc(), name());
282 log << MSG::INFO << "in excute()" << endreq;
283
284
285 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
286 if (!eventHeader) {
287 log << MSG::FATAL << "Could not find Event Header" << endreq;
288 return StatusCode::FAILURE;
289 }
290 evt = eventHeader->eventNumber();
291
292 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
293 if (!recCgemSegmentCol){
294 log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
295 unSegNum++;
296 return StatusCode::SUCCESS;
297 }
298
299 /// Get RawDataProviderSvc
300 IRawDataProviderSvc* irawDataProviderSvc;
301 StatusCode sc = service ("RawDataProviderSvc", irawDataProviderSvc);
302 RawDataProviderSvc *_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
303 if ( sc.isFailure() ){
304 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
305 return StatusCode::FAILURE;
306 }
307
308 uint32_t getDigiFlag = 0;
309 if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
310 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
311 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
312 MdcDigiVec mdcDigiVec = _rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
313 if (0 == mdcDigiVec.size()){
314 log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
315 // return StatusCode::SUCCESS;
316 }
317
318 MdcDigiVec::iterator itDigi = mdcDigiVec.begin();
319
320
321 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
322 if (!recMdcTrackCol){
323 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
324 unTrkNum++;
325 return StatusCode::SUCCESS;
326 }
327
328 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
329 if (!mcParticleCol) {
330 log << MSG::ERROR<< "Could not find McParticle" << endreq;
331 } else{
332 if(m_debug==1) cout<< "size_mcParticleCol "<<mcParticleCol->size()<< endl;
333 McParticleCol::iterator iter_mc = mcParticleCol->begin();
334 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
335 bool findJpsi=false;
336 int pdg=(*iter_mc)->particleProperty();
337 if(abs(pdg)==13){
338 findJpsi=true;
339 }
340 if(!findJpsi) continue;
341 double mc_charge;
342
343 if(pdg<0)mc_charge = 1.;
344 else mc_charge =-1.;
345 double truth_x = (*iter_mc)->initialPosition().x();
346 double truth_y = (*iter_mc)->initialPosition().y();
347 double truth_z = (*iter_mc)->initialPosition().z();
348
349 double px = (*iter_mc)->initialFourMomentum().px();
350 double py = (*iter_mc)->initialFourMomentum().py();
351 double pz = (*iter_mc)->initialFourMomentum().pz();
352
353 HepPoint3D truth_position(truth_x, truth_y, truth_z);
354 Hep3Vector p3(px, py, pz);
355 KalmanFit::Helix* tmp_helix = new KalmanFit::Helix(truth_position, p3, mc_charge);
356 HepPoint3D tmp_pivot(0., 0., 0.);
357 tmp_helix->pivot(tmp_pivot);
358 truth_dr = tmp_helix->dr();
359 truth_dz = tmp_helix->dz();
360 truth_phi0 = tmp_helix->phi0();
361 truth_kappa = tmp_helix->kappa();
362 truth_tanl = tmp_helix->tanl();
363
364 costhetaMC = pz/sqrt(pow(px,2)+pow(py,2)+pow(pz,2));
365 phiMC = p3.phi();
366 pdgMC = pdg;
367
368 delete tmp_helix;
369 }
370 }
371
372
373
374 for(;itDigi != mdcDigiVec.end(); itDigi++){
375 Identifier mdcId = (*itDigi)->identify();
376 long l = MdcID::layer(mdcId);
377 long w = MdcID::wire(mdcId);
378 }
379
380 int trk = 0;
381 RecCgemSegmentCol::iterator itSeg = recCgemSegmentCol->begin();
382 for(match_nsegment = 0, segment_size = 0;itSeg != recCgemSegmentCol->end(); itSeg++, segment_size++){
383 RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
384 if((*itSeg)->gethelix(1) < -900) match_failed[segment_size] = 0;
385 else match_failed[segment_size] = 1;
386 for(track_size = 0;itTrk != recMdcTrackCol->end(); itTrk++, track_size++, match_nsegment++){
387 //if(!track_size){cout<<"Track size = "<<track_size<<endl;trk++;}
388 match_dr_chisq[match_nsegment] = GetSigChisq((*itSeg)->gethelix(0), (*itTrk)->helix(0), (*itSeg)->gethelix_err(0) + (*itTrk)->err(0));
389 match_phi0_chisq[match_nsegment] = GetSigChisq((*itSeg)->gethelix(1), (*itTrk)->helix(1), (*itSeg)->gethelix_err(5) + (*itTrk)->err(5));
390 match_kappa_chisq[match_nsegment] = GetSigChisq((*itSeg)->gethelix(2), (*itTrk)->helix(2), (*itSeg)->gethelix_err(9) + (*itTrk)->err(9));
391 match_dz_chisq[match_nsegment] = GetSigChisq((*itSeg)->gethelix(3), (*itTrk)->helix(3), (*itSeg)->gethelix_err(12) + (*itTrk)->err(12));
392 match_tanl_chisq[match_nsegment] = GetSigChisq((*itSeg)->gethelix(4), (*itTrk)->helix(4), (*itSeg)->gethelix_err(14) + (*itTrk)->err(14));
393
394 double cut[5];double pt_tmp;
395 if((*itTrk)->helix(2) != 0) pt_tmp = 1./(*itTrk)->helix(2);
396 else{
397 log << MSG::FATAL << "The p_t of track is 0!!!" << endreq;
398 return StatusCode::FAILURE;
399 }
400 cut[0] = GetWindowChisq(pt_tmp, m_Pt, m_w_dr);
401 cut[1] = GetWindowChisq(pt_tmp, m_Pt, m_w_phi0);
402 cut[2] = GetWindowChisq(pt_tmp, m_Pt, m_w_kappa);
403 cut[3] = GetWindowChisq(pt_tmp, m_Pt, m_w_dz);
404 cut[4] = GetWindowChisq(pt_tmp, m_Pt, m_w_tanl);
405
406 // pt_tmp = fabs((*itTrk)->helix(2));
407 // cut[0] = GetWindowChisq(pt_tmp, m_dr_par[0], m_dr_par[1], m_dr_par[2], m_dr_par[3]);
408 // cut[1] = GetWindowChisq(pt_tmp, m_phi0_par[0], m_phi0_par[1], m_phi0_par[2], m_phi0_par[3]);
409 // cut[2] = GetWindowChisq(pt_tmp, m_kappa_par[0], m_kappa_par[1], m_kappa_par[2], m_kappa_par[3]);
410 // cut[3] = GetWindowChisq(pt_tmp, m_dz_par[0], m_dz_par[1], m_dz_par[2], m_dz_par[3]);
411 // cut[4] = GetWindowChisq(pt_tmp, m_tanl_par[0], m_tanl_par[1], m_tanl_par[2], m_tanl_par[3]);
412 // cout<<"chisq_dr = "<<match_dr_chisq[match_nsegment]<<"\t"<<cut[0]<<endl<<
413 // "chisq_phi0 = "<<match_phi0_chisq[match_nsegment]<<"\t"<<cut[1]<<endl<<
414 // "chisq_kappa = "<<match_kappa_chisq[match_nsegment]<<"\t"<<cut[2]<<endl<<
415 // "chisq_dz = "<<match_dz_chisq[match_nsegment]<<"\t"<<cut[3]<<endl<<
416 // "chisq_tanl = "<<match_tanl_chisq[match_nsegment]<<"\t"<<cut[4]<<endl<<endl;
417 if(match_dr_chisq[match_nsegment]>cut[0]||
418 match_phi0_chisq[match_nsegment]>cut[1]||
419 match_kappa_chisq[match_nsegment]>cut[2]||
420 match_dz_chisq[match_nsegment]>cut[3]||
421 match_tanl_chisq[match_nsegment]>cut[4])EvOutWindow++;
422 else EvInWindow++;
423
424 match_IsTruth[match_nsegment] = (*itSeg)->getmatch();
425
426 track_dr[match_nsegment] = (*itTrk)->helix(0);
427 track_phi0[match_nsegment] = (*itTrk)->helix(1);
428 track_kappa[match_nsegment] = (*itTrk)->helix(2);
429 track_dz[match_nsegment] = (*itTrk)->helix(3);
430 track_tanl[match_nsegment] = (*itTrk)->helix(4);
431 track_dr_err[match_nsegment] = (*itTrk)->err(0);
432 track_phi0_err[match_nsegment] = (*itTrk)->err(5);
433 track_kappa_err[match_nsegment] = (*itTrk)->err(9);
434 track_dz_err[match_nsegment] = (*itTrk)->err(12);
435 track_tanl_err[match_nsegment] = (*itTrk)->err(14);
436 seg_dr[match_nsegment] = (*itSeg)->gethelix(0);
437 seg_phi0[match_nsegment] = (*itSeg)->gethelix(1);
438 seg_kappa[match_nsegment] = (*itSeg)->gethelix(2);
439 seg_dz[match_nsegment] = (*itSeg)->gethelix(3);
440 seg_tanl[match_nsegment] = (*itSeg)->gethelix(4);
441 seg_dr_err[match_nsegment] = (*itSeg)->gethelix_err(0);
442 seg_phi0_err[match_nsegment] = (*itSeg)->gethelix_err(5);
443 seg_kappa_err[match_nsegment] = (*itSeg)->gethelix_err(9);
444 seg_dz_err[match_nsegment] = (*itSeg)->gethelix_err(12);
445 seg_tanl_err[match_nsegment] = (*itSeg)->gethelix_err(14);
446
447 //int trackID = (*itTrk)->trackId();
448 //HitRefVec vechits = (*itTrk)->getVecHits();
449 //for(int i = 0; i<vechits.size(); i++){
450 // Identifier mdcId = (vechits[i])->getMdcId();
451 // long l = MdcID::layer(mdcId);
452 // long w = MdcID::wire(mdcId);
453 //}
454 }
455 }
456 //cout<<" trk = " << trk<<endl;
457
458 match_ntuple->write();
459
460 return StatusCode::SUCCESS;
461
462}
463
464StatusCode CgemMdcCombAlg::exe_v2()
465{
466 MsgStream log(msgSvc(), name());
467 log << MSG::INFO << "in excute()" << endreq;
468 if(m_debug) cout<<"--- CgemMdcCombAlg::exe ---"<<endl;
469
470 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
471 if (!eventHeader) {
472 log << MSG::FATAL << "Could not find Event Header" << endreq;
473 return StatusCode::FAILURE;
474 }
475 m_run = eventHeader->runNumber();
476 if(m_run<0&&m_prodNT) getMcPar();
477 m_evt = eventHeader->eventNumber();
478
479 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
480 if (!recCgemClusterCol){
481 log << MSG::WARNING << "Could not retrieve Cgem Cluster list" << endreq;
482 }
483
484 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
485 int nCgemSeg = 0;
486 if (recCgemSegmentCol){
487 nCgemSeg = recCgemSegmentCol->size();
488 }
489 else {
490 log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
491 return StatusCode::SUCCESS;
492 }
493 if(m_debug) cout<<"nCgemSeg = "<<nCgemSeg<<endl;
494 //if(nCgemSeg==0) return StatusCode::SUCCESS;
495 RecCgemSegmentCol::iterator itCgemSeg = recCgemSegmentCol->begin();
496
497 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
498 int nMDCTrk = 0;
499 if (recMdcTrackCol){
500 nMDCTrk = recMdcTrackCol->size();
501 }
502 else {
503 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
504 return StatusCode::SUCCESS;
505 }
506 if(m_debug) cout<<"nMDCTrk = "<<nMDCTrk<<endl;
507 //if(nMDCTrk==0) return StatusCode::SUCCESS;
508
509 //if(m_debug) {
510 // cout<<"Before match:"<<endl;
511 // printMdcTrk();
512 //}
513
514 // temp
515 //double sigma_par[5] = {0.103, 0.0161, 0.4246, 0.502, 0.011};
516 //double sigma_par[5] = {0.103, 0.0161, 0.4246, 2, 0.05};
517
518 if(m_prodNT) {match_nsegment=0;}
519 if(nMDCTrk>0)
520 {
521 RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
522 int i_ODC = 0;
523 for(;itTrk != recMdcTrackCol->end(); itTrk++, i_ODC++)// loop ODC tracks
524 {
525 RecCgemSegmentCol::iterator itCgemSegMatched;
526 RecCgemSegmentCol::iterator itCgemSeg = recCgemSegmentCol->begin();
527 for(; itCgemSeg!=recCgemSegmentCol->end(); itCgemSeg++)// loop seg
528 {
529 int i_err=0;
530 double pt=0.5;
531 if((*itTrk)->helix(2) != 0) pt=1./(*itTrk)->helix(2);
532 bool match=true;
533 bool matched[5]={true, true, true, true, true};
534 bool matchSig = true;
535 bool matchedSig[5]={true, true, true, true, true};
536 double totChi2(0), totChi2InXY(0), totChi2PhiTheta(0), totChi2PhiThetaKappa(0);
537 double chi2_Sig[5] = {9999, 9999, 9999, 9999, 9999};
538 for(int i=0; i<5; i++)
539 {
540 //cout<<"i_err = "<<i_err<<endl;
541 double parSeg = (*itCgemSeg)->gethelix(i);
542 double parTrk = (*itTrk)->helix(i);
543 //if(fabs(parSeg-parTrk)>5.0*sigma_par[i])
544 double delta=parSeg-parTrk;
545 if(i==1) {
546 while(delta> CLHEP::pi) delta-=CLHEP::twopi;
547 while(delta<-CLHEP::pi) delta+=CLHEP::twopi;
548 }
549 double chi2_new = Chi2Match(delta, pt, i);// chi2 after calibration of sigma
550 chi2_Sig[i]=chi2_new;
551 if(chi2_new>mv_matchChi2Cut[i])
552 {
553 matchSig=false;
554 matchedSig[i]=false;
555 }
556 if(i==1||i==4||i==2) totChi2PhiThetaKappa+=chi2_new;
557 double errSqTot = (*itCgemSeg)->gethelix_err(i_err)+(*itTrk)->err(i_err);
558 double chi2Match = GetSigChisq(parSeg, parTrk, errSqTot);
559 totChi2+=chi2Match;
560 if(i<3) totChi2InXY+=chi2Match;
561 if(i==1||i==4) totChi2PhiTheta+=chi2Match;
562 double cut = GetWindowChisq(pt, i);
563 if(chi2Match>cut)
564 {
565 match=false;
566 matched[i]=false;
567 }
568 if(m_debug) cout<<"i, seg, trk, chi2Match = "<<i<<", "<<parSeg<<", "<<parTrk<<", "<<chi2Match<<endl;
569 i_err += 5-i;
570 if(m_prodNT&&match_nsegment<1000) {
571 switch(i) {
572 case 0:
573 match_dr_chisq[match_nsegment]=chi2Match;
574 track_dr[match_nsegment]=parTrk;
575 seg_dr[match_nsegment]=parSeg;
576 break;
577 case 1:
578 match_phi0_chisq[match_nsegment]=chi2Match;
579 track_phi0[match_nsegment]=parTrk;
580 seg_phi0[match_nsegment]=parSeg;
581 break;
582 case 2:
583 match_kappa_chisq[match_nsegment]=chi2Match;
584 track_kappa[match_nsegment]=parTrk;
585 seg_kappa[match_nsegment]=parSeg;
586 break;
587 case 3:
588 match_dz_chisq[match_nsegment]=chi2Match;
589 track_dz[match_nsegment]=parTrk;
590 seg_dz[match_nsegment]=parSeg;
591 break;
592 case 4:
593 match_tanl_chisq[match_nsegment]=chi2Match;
594 track_tanl[match_nsegment]=parTrk;
595 seg_tanl[match_nsegment]=parSeg;
596 break;
597 default:
598 break;
599 }
600 }
601 }// loop Helix 5 parameters
602 if(m_prodNT&&match_nsegment<1000) match_nsegment++;
603 bool matchedInXY = matched[0]*matched[1]*matched[2];
604 bool matchPhiTheta = matched[1]*matched[4];
605 if(m_debug) {
606 cout<<"matched "<<match<<endl;
607 cout<<"matched only in r-phi: "<<matchedInXY<<endl;
608 cout<<"matched only in phi-theta: "<<matchPhiTheta<<endl;
609 cout<<"match with sigma "<<matchSig<<endl;
610 cout<<"match with sigma [] "<<matchedSig[0]<<matchedSig[1]<<matchedSig[2]<<matchedSig[3]<<matchedSig[4]<<endl;
611 cout<<"chi2 match with sigma [] "<<chi2_Sig[0]<<", "<<chi2_Sig[1]<<", "<<chi2_Sig[2]<<", "<<chi2_Sig[3]<<", "<<chi2_Sig[4]<<endl;
612 }
613 bool matchFlag = matchSig;// options: match, matchSig,
614 double matchChi2 = totChi2PhiThetaKappa;// options: totChi2PhiTheta, totChi2PhiThetaKappa
615 if(matchFlag)
616 {
617 if(m_prodNT) CgemSegMatched[i_ODC]=1;
618 //cout<<"matchChi2: "<<(*itTrk)->genMatchChi2()<<endl;
619 if( (*itTrk)->getNcluster() == 0 ||
620 ((*itTrk)->getNcluster()!= 0 && (*itTrk)->getMatchChi2()>matchChi2) )
621 {
622 ClusterRefVec vecclusters;
623 RecCgemClusterCol::iterator itClusterBegin = recCgemClusterCol->begin();
624 int trkid = (*itTrk)->trackId();
625 (*(itClusterBegin+(*itCgemSeg)->getclusterid_1()))->setTrkId(trkid);
626 (*(itClusterBegin+(*itCgemSeg)->getclusterid_2()))->setTrkId(trkid);
627 (*(itClusterBegin+(*itCgemSeg)->getclusterid_3()))->setTrkId(trkid);
628 vecclusters.push_back(*(itClusterBegin+(*itCgemSeg)->getclusterid_1()));
629 vecclusters.push_back(*(itClusterBegin+(*itCgemSeg)->getclusterid_2()));
630 vecclusters.push_back(*(itClusterBegin+(*itCgemSeg)->getclusterid_3()));
631 (*itTrk)->setVecClusters(vecclusters);
632 (*itTrk)->setNcluster(vecclusters.size());
633 (*itTrk)->setMatchChi2(matchChi2);
634 (*itTrk)->setNdof( (*itTrk)->ndof() + vecclusters.size() );
635 (*itTrk)->setNster( (*itTrk)->nster() + vecclusters.size() );
636 (*itTrk)->setFirstLayer(0);
637 (*itTrk)->setNlayer( (*itTrk)->nlayer() + vecclusters.size() );
638 (*itTrk)->setNCgemXCluster(vecclusters.size());
639 (*itTrk)->setNCgemVCluster(vecclusters.size());
640 itCgemSegMatched = itCgemSeg;
641 }
642 }
643 //else if(m_prodNT) CgemSegMatched[i_ODC]=0;
644 }// loop seg
645 // average dr dz tanL --> begin
646 if((*itTrk)->getNcluster()!= 0)
647 {
648 HepVector a_trk = (*itTrk)->helix();
649 int iPar[3]={0, 3, 4};
650 HepSymMatrix parTrkErMat = (*itTrk)->err();
651 //cout<<"Track error matrix: "<<parTrkErMat<<endl;
652 for(int i=0; i<3; i++)
653 {
654 int j = iPar[i];
655 double parSeg = (*itCgemSegMatched)->gethelix(j);
656 double parSegEr = (*itCgemSegMatched)->gethelix_err(j*(11-j)/2);
657 double parTrk = (*itTrk)->helix(j);
658 double parTrkEr = parTrkErMat(j+1,j+1); //(*itTrk)->err(j);
659 //cout<<"i, errCGEM, errMDC = "<<j<<", "<<parSegEr<<", "<<parTrkEr<<endl;
660 double parAverage = (parSeg/parSegEr+parTrk/parTrkEr)/(1./parSegEr+1./parTrkEr);
661 //cout<<"parSeg, parTrk, parAverage = "<<parSeg<<", "<<parTrk<<", "<<parAverage<<endl;
662 a_trk[j]=parAverage;
663 if(m_prodNT)
664 {
665 if(i==0) dr_avg[i_ODC]=parAverage;
666 if(i==1) dz_avg[i_ODC]=parAverage;
667 if(i==2) tanL_avg[i_ODC]=parAverage;
668 }
669 }
670 (*itTrk)->setHelix(a_trk);
671 }
672 else {
673 if(m_prodNT)
674 {
675 dr_avg[i_ODC]=(*itTrk)->helix(0);
676 dz_avg[i_ODC]=(*itTrk)->helix(3);
677 tanL_avg[i_ODC]=(*itTrk)->helix(4);
678 }
679 }
680 // average dr dz tanL --> end
681
682 }// loop track
683 }// if nMDC>0
684
685 if(m_debug) {
686 cout<<"After match:"<<endl;
687 printMdcTrk();
688 }
689
690 if(m_prodNT) {
691 run = m_run;
692 evt = m_evt;
693 track_size=nMDCTrk;
694 segment_size=nCgemSeg;
695 match_ntuple->write();
696 }
697
698 return StatusCode::SUCCESS;
699}
700
701double CgemMdcCombAlg::GetSigChisq(double x, double y, double err) {
702 return pow(x -y, 2)/err;
703}
704
705double CgemMdcCombAlg::GetWindowChisq(double pt, double a0, double a1, double a2, double a3) {
706 return a0+a1*pt+a2*pt*pt+a3*pt*pt*pt;
707}
708
709double CgemMdcCombAlg::GetWindowChisq(double pt, vector<double> Pt, vector<double> w_p) {
710 int n = Pt.size();
711 double tmp_pt[n]; double tmp_p[n];
712 //int n = 7;
713 for(int mm = 0; mm!=n; mm++){
714 tmp_pt[mm] = Pt[mm];
715 tmp_p[mm] = w_p[mm];
716 }
717 TGraph *tmp_w = new TGraph(n,tmp_pt,tmp_p);
718 double y = tmp_w->Eval(fabs(pt));
719 //cout<<"pt = " << pt<<"\ty = "<<y<<endl;
720 delete tmp_w;
721 return y;
722}
723
724double CgemMdcCombAlg::GetWindowChisq(double pt, int i_par)
725{
726 vector<double> w_par;
727 switch(i_par)
728 {
729 case 0:
730 w_par = m_w_dr;
731 break;
732 case 1:
733 w_par = m_w_phi0;
734 break;
735 case 2:
736 w_par = m_w_kappa;
737 break;
738 case 3:
739 w_par = m_w_dz;
740 break;
741 case 4:
742 w_par = m_w_tanl;
743 break;
744
745 default:
746 cout<<"CgemMdcCombAlg::GetWindowChisq() wrong i_par! "<<endl;
747 return 9999.;
748 break;
749 }
750
751 return GetWindowChisq(pt, m_Pt, w_par);
752
753}
754
755
756double CgemMdcCombAlg::GetTotChisq(double seg[5], double trk[5], double err[15]) {
757 HepMatrix delta_Par(5,1);
758 HepSymMatrix ParErr(5);
759 int m = 0;
760 for(int i = 0;i<5;i++){
761 delta_Par(i+1,1) = seg[i] - trk[i];
762 for(int j = i; j<5;j++){
763 ParErr(i+1,j+1) = err[m];
764 if(i!=j)ParErr(j+1,i+1) = ParErr(i+1,j+1);
765 m++;
766 }
767 }
768 int ierr = 1;
769 ParErr.invert(ierr);
770 double chisq = (ParErr.similarityT(delta_Par)).determinant();
771 if(!ierr) return chisq;
772 else return -999.;
773}
774
775double CgemMdcCombAlg::Chi2Match(double delta, double pt, int iPar)
776{
777 pt=fabs(pt);
778 //if(pt<mv_Pt[0]) pt=mv_Pt[0];
779 if(pt>mv_Pt[m_nPt-1]) pt=mv_Pt[m_nPt-1];
780
781 TGraph *g_mean;
782 TGraph *g_sig;
783 switch(iPar)
784 {
785 case 0:
786 g_mean = mg_mean_dr;
787 g_sig = mg_sig_dr;
788 break;
789 case 1:
790 g_mean = mg_mean_phi0;
791 g_sig = mg_sig_phi0;
792 break;
793 case 2:
794 g_mean = mg_mean_kappa;
795 g_sig = mg_sig_kappa;
796 break;
797 case 3:
798 g_mean = mg_mean_dz;
799 g_sig = mg_sig_dz;
800 break;
801 case 4:
802 g_mean = mg_mean_tanL;
803 g_sig = mg_sig_tanL;
804 break;
805
806 default:
807 cout<<"CgemMdcCombAlg::Chi2Match() wrong iPar! "<<endl;
808 return 9999.;
809 break;
810 }
811 double mean = g_mean->Eval(pt);
812 double sigma = g_sig->Eval(pt);
813 double chi2 = pow((delta-mean)/sigma,2.0);
814 //cout<<"delta, pt, iPar, = "<<delta<<", "<<pt<<", "<<iPar<<", "<<endl;
815 //cout<<"mean, sigma, chi2 = "<<mean<<", "<<sigma<<", "<<chi2<<endl;
816 return chi2;
817}
818
819// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
820StatusCode CgemMdcCombAlg::finalize() {
821
822 MsgStream log(msgSvc(), name());
823 log << MSG::INFO << "in finalize()" << endreq;
824
825 /*
826 cout<<"segment leave number = "<<unSegNum<<endl;
827 cout<< "track leave number = "<<unTrkNum<<endl;
828 cout<<"events in window = " <<EvInWindow<<endl;
829 cout<<"events out window = " <<EvOutWindow<<endl;
830 */
831 if(m_nPt) delete mg_mean_dr;
832 if(mg_mean_phi0) delete mg_mean_phi0;
833 delete mg_mean_kappa;
834 delete mg_mean_dz;
835 delete mg_mean_tanL;
836 delete mg_sig_dr;
837 delete mg_sig_phi0;
838 delete mg_sig_kappa;
839 delete mg_sig_dz;
840 delete mg_sig_tanL;
841
842 return StatusCode::SUCCESS;
843}
844
845void CgemMdcCombAlg::printMdcTrk()
846{
847 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
848 //SmartDataPtr<DstMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Dst/DstMdcTrackCol");
849 if (!recMdcTrackCol){
850 cout<<"printMdcTrk() can not find RecMdcTrackCol"<<endl;
851 return;
852 }
853 else {
854 int nMDCTrk = recMdcTrackCol->size();
855 RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
856 //DstMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
857 int i_ODC = 0;
858 for(;itTrk != recMdcTrackCol->end(); itTrk++, i_ODC++)// loop ODC tracks
859 {
860 int ndof = (*itTrk)->ndof();
861 int nster = (*itTrk)->nster();
862 int nlayer = (*itTrk)->nlayer();
863 int firstLayerId = (*itTrk)->firstLayer();
864 int lastLayerId = (*itTrk)->lastLayer();
865 int nClu2D = (*itTrk)->getNcluster();
866 int nCgemCluX = (*itTrk)->nCgemXCluster();
867 int nCgemCluV = (*itTrk)->nCgemVCluster();
868 if(i_ODC==0)
869 cout <<setw(15)<<"ndof"<<setw(15)<<"nster"<<setw(15)<<"nlayer"<<setw(15)<<"firstLayerId"<<setw(15)<<"lastLayerId"
870 <<setw(15)<<"nClu2D"
871 <<setw(15)<<"nCgemCluX"<<setw(15)<<"nCgemCluV"<<endl;
872 cout<<setw(15)<<ndof<<setw(15)<<nster<<setw(15)<<nlayer<<setw(15)<<firstLayerId<<setw(15)<<lastLayerId
873 <<setw(15)<<nClu2D
874 <<setw(15)<<nCgemCluX<<setw(15)<<nCgemCluV<<endl;
875 }
876 }
877}
n
const Int_t n
Definition: DataBase/tau_mode.c:65
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
cut
*********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 cut
Definition: KarFin.h:27
MdcDigiVec
std::vector< MdcDigi * > MdcDigiVec
Definition: MdcRawDataProvider.h:11
ClusterRefVec
SmartRefVector< RecCgemCluster > ClusterRefVec
Definition: RecMdcTrack.h:28
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
CgemMdcCombAlg::CgemMdcCombAlg
CgemMdcCombAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition: CgemMdcCombAlg.cxx:67
CgemMdcCombAlg::initialize
StatusCode initialize()
Definition: CgemMdcCombAlg.cxx:104
CgemMdcCombAlg::finalize
StatusCode finalize()
Definition: CgemMdcCombAlg.cxx:820
CgemMdcCombAlg::execute
StatusCode execute()
Definition: CgemMdcCombAlg.cxx:223
KalmanFit::Helix::dr
double dr(void) const
returns an element of parameters.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:229
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:211
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
RawDataProviderSvc::getMdcDigiVec
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
Definition: RawDataProviderSvc.h:29
Event
Definition: Event.h:21