CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemSegmentFitAlg.cxx
Go to the documentation of this file.
1#include "CgemSegmentFitAlg/CgemSegmentFitAlg.h"
2#include "CgemSegmentFitAlg/CgemSegmentFun.h"
3#include "GaudiKernel/IDataManagerSvc.h"
4#include "GaudiKernel/SmartDataPtr.h"
5#include "GaudiKernel/IDataProviderSvc.h"
6#include "GaudiKernel/IPartPropSvc.h"
7#include "GaudiKernel/MsgStream.h"
8#include "GaudiKernel/AlgFactory.h"
9#include "GaudiKernel/ISvcLocator.h"
10#include "GaudiKernel/DataSvc.h"
11#include "GaudiKernel/PropertyMgr.h"
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/INTupleSvc.h"
14
15
16#include "EventModel/EventHeader.h"
17#include "CLHEP/Units/PhysicalConstants.h"
18
19#include "HepPDT/ParticleDataTable.hh"
20#include "HepPDT/ParticleData.hh"
21
22#include "EventModel/Event.h"
23#include "EvTimeEvent/RecEsTime.h"
24#include "EvtRecEvent/EvtRecEvent.h"
25#include "EvtRecEvent/EvtRecTrack.h"
26#include "Identifier/Identifier.h"
27#include "CLHEP/Vector/ThreeVector.h"
28#include "CLHEP/Matrix/Vector.h"
29#include "CLHEP/Matrix/SymMatrix.h"
30#include "CLHEP/Geometry/Point3D.h"
31#include "CLHEP/Vector/LorentzVector.h"
32//#include "TrkFitter/TrkHelixMaker.h"
33#include "KalFitAlg/helix/Helix.h"
34
35#include "CLHEP/Matrix/Matrix.h"
36
37#include "CgemGeomSvc/CgemGeomSvc.h"
38#include "CgemRecEvent/RecCgemCluster.h"
39#include "CgemRecEvent/RecCgemSegment.h"
40#include "CgemClusterCreate/CgemClusterCreate.h"
41#include "McTruth/CgemMcHit.h"
42#include "CgemRawEvent/CgemDigi.h"
43
44
45#include "TRandom.h"
46#include "TArrayI.h"
47#include "TGraph.h"
48#include "TF1.h"
49#include "TMath.h"
50#include "TMinuit.h"
51
52
53#include <vector>
54#include <iostream>
55#include <fstream>
56#include <iomanip>
57#include <cstring>
58
59using namespace std;
60using namespace Event;
61//using namespace KalmanFit;
62
63/////////////////////////////////////////////////////////////////////////////
64CgemSegmentFitAlg::CgemSegmentFitAlg(const std::string& name, ISvcLocator* pSvcLocator) :
65 Algorithm(name, pSvcLocator)
66{
67 declareProperty("debug", m_debug = 0);
68 declareProperty("check", m_checkIdx = 0);
69 declareProperty("version", m_ver = 2);
70}
71
72// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
73StatusCode CgemSegmentFitAlg::initialize(){
74
75 MsgStream log(msgSvc(), name());
76 log << MSG::INFO << "in initialize()" << endreq;
77 if(m_checkIdx>0) {
78 /////////////////////////////// book ntuple //////////////////////////////////////
79 INTupleSvc* ntupleSvc;
80 StatusCode sc = service("NTupleSvc", ntupleSvc);
81 NTuplePtr mt(ntupleSvc, "FILE169/CgemSegmentHelix");
82 if ( mt ) helix_ntuple = mt;
83 else{
84 helix_ntuple = ntupleSvc->book("FILE169/CgemSegmentHelix", CLID_ColumnWiseTuple, "SegFit");
85 StatusCode status;
86 if(helix_ntuple){
87 status = helix_ntuple->addItem("run", cgem_run);
88 status = helix_ntuple->addItem("evt", cgem_evt);
89 status = helix_ntuple->addItem("ntrack",truth_ntrack,0,10);
90 status = helix_ntuple->addIndexedItem("truth_dr", truth_ntrack , truth_dr);
91 status = helix_ntuple->addIndexedItem("truth_dz", truth_ntrack , truth_dz);
92 status = helix_ntuple->addIndexedItem("truth_phi0", truth_ntrack , truth_phi0);
93 status = helix_ntuple->addIndexedItem("truth_kappa", truth_ntrack , truth_kappa);
94 status = helix_ntuple->addIndexedItem("truth_tanl", truth_ntrack , truth_tanl);
95 status = helix_ntuple->addIndexedItem("truth_charge", truth_ntrack , truth_charge);
96 status = helix_ntuple->addIndexedItem("truth_CosTheta", truth_ntrack , truth_CosTheta);
97 status = helix_ntuple->addItem("nsegment",cgem_nsegment,0,1000);
98 status = helix_ntuple->addIndexedItem("segmentid", cgem_nsegment , cgem_segmentid);
99 status = helix_ntuple->addIndexedItem("dr_ini", cgem_nsegment , cgem_dr_ini);
100 status = helix_ntuple->addIndexedItem("dz_ini", cgem_nsegment , cgem_dz_ini);
101 status = helix_ntuple->addIndexedItem("phi0_ini", cgem_nsegment , cgem_phi0_ini);
102 status = helix_ntuple->addIndexedItem("kappa_ini", cgem_nsegment , cgem_kappa_ini);
103 status = helix_ntuple->addIndexedItem("tanl_ini", cgem_nsegment , cgem_tanl_ini);
104 status = helix_ntuple->addIndexedItem("dr", cgem_nsegment , cgem_dr);
105 status = helix_ntuple->addIndexedItem("dz", cgem_nsegment , cgem_dz);
106 status = helix_ntuple->addIndexedItem("phi0", cgem_nsegment , cgem_phi0);
107 status = helix_ntuple->addIndexedItem("kappa", cgem_nsegment , cgem_kappa);
108 status = helix_ntuple->addIndexedItem("tanl", cgem_nsegment , cgem_tanl);
109 status = helix_ntuple->addIndexedItem("chisq", cgem_nsegment , cgem_chisq);
110 status = helix_ntuple->addIndexedItem("fitFlag", cgem_nsegment , cgem_fitFlag);
111 status = helix_ntuple->addIndexedItem("charge", cgem_nsegment , cgem_charge);
112
113 status = helix_ntuple->addItem("n_layer",n_layer,0,3);
114 status = helix_ntuple->addIndexedItem("n_stripX", n_layer , n_stripX);
115 status = helix_ntuple->addIndexedItem("n_stripV", n_layer , n_stripV);
116 status = helix_ntuple->addIndexedItem("mc_phi", n_layer , mc_phi);
117 status = helix_ntuple->addIndexedItem("mc_v", n_layer , mc_v);
118 status = helix_ntuple->addIndexedItem("rec_phi", n_layer , rec_phi);
119 status = helix_ntuple->addIndexedItem("rec_v", n_layer , rec_v);
120 }
121 }
122 }
123 check_segNum = 0;
124 check_segOver = 0;
125 wide_delta = 570./2/10000; //cm
126 pitch = (650.-570.)/2/10000; //cm
127
128 ISvcLocator* svcLocator = Gaudi::svcLocator();
129 ICgemGeomSvc* ISvc;
130 StatusCode sc=svcLocator->service("CgemGeomSvc", ISvc);
131 CgemGeomSvc* myCgemGeomSvc=dynamic_cast<CgemGeomSvc *>(ISvc);
132 int nCgemLayer = myCgemGeomSvc->getNumberOfCgemLayer();
133 for(int i=0; i<nCgemLayer; i++)
134 {
135 CgemGeoLayer* CgemLayer = myCgemGeomSvc->getCgemLayer(i);
136 CgemSegmentFun::rl[i]=(CgemLayer->getMiddleROfGapD())/10.0;//cm
137 CgemSegmentFun::a_stero[i]=(CgemLayer->getAngleOfStereo())/180.0*CLHEP::pi;// degree -> radians
138 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(i,0);
139 CgemSegmentFun::r_X[i]=readoutPlane->getRX()/10;//cm
140 CgemSegmentFun::r_V[i]=readoutPlane->getRV()/10;//cm
141 if(m_debug) {
142 cout<<"iLayer, r, ang_stereo, Rx, Rv = "<<i<<", "<<CgemSegmentFun::rl[i]<<", "<<CgemSegmentFun::a_stero[i]<<", "<<CgemSegmentFun::r_X[i]<<", "<<CgemSegmentFun::r_V[i]<<endl;
143 }
144 }
145
146 sc = service ("CgemCalibFunSvc", myCgemCalibSvc);
147 if ( sc.isFailure() ){
148 cout<< name() << "Could not load MdcCalibFunSvc!" << endl;
149 return sc;
150 }
151
152 return StatusCode::SUCCESS;
153
154}
155
156StatusCode CgemSegmentFitAlg::execute()
157{
158 if(m_ver==1) exe_v1();
159 else if(m_ver==2) exe_v2();
160
161 return StatusCode::SUCCESS;
162}
163
164StatusCode CgemSegmentFitAlg::exe_v1()
165{
166 MsgStream log(msgSvc(), name());
167 log << MSG::INFO << "in excute()" << endreq;
168
169 const double pi = 4.*atan(1);
170 n_layer = 3;
171
172 // SmartDataPtr<CgemGeomSvc> cgemGeom(eventSvc(), "CgemGeomSvc");
173 // if (!cgemGeom){
174 // log << MSG::WARNING << "Could not retrieve Cgem Geomtry list" << endreq;
175 // return StatusCode::FAILURE;
176 // }
177 // CgemGeomSvc::iterator itGeom = cgemGeom->begin();
178
179 int vec_clusterid[3];
180 m_seg_map.clear();
181
182 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
183 if (!eventHeader) {
184 log << MSG::FATAL << "Could not find Event Header" << endreq;
185 return StatusCode::FAILURE;
186 }
187 cgem_evt = eventHeader->eventNumber();
188 cgem_run = eventHeader->runNumber();
189
190 // if(cgem_evt<=221) return StatusCode::FAILURE;
191 // cout<<"event_ID = "<<cgem_evt<<endl;
192 //if(cgem_evt%1000==0)cout<<"event_ID = "<<cgem_evt<<endl;
193 // if(cgem_evt > 62)return StatusCode::SUCCESS;
194
195 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
196 if (!recCgemSegmentCol){
197 log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
198 check_segNum++;
199 // return StatusCode::FAILURE;
200 return StatusCode::SUCCESS;
201 }
202 RecCgemSegmentCol::iterator itSeg;
203 /*
204 itSeg = recCgemSegmentCol->begin();
205 for(;itSeg != recCgemSegmentCol->end(); itSeg++){
206 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
207 cout<<"event ID = "<<(*itSeg)->geteventid()<<endl;;
208 cout<<"segment ID = "<<(*itSeg)->getsegmentid()<<endl;;
209 cout<<"cluster ID_1 = "<<(*itSeg)->getclusterid_1()<<endl;;
210 cout<<"cluster ID_2 = "<<(*itSeg)->getclusterid_2()<<endl;;
211 cout<<"cluster ID_3 = "<<(*itSeg)->getclusterid_3()<<endl;;
212 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
213 }
214 */
215 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
216 if (!recCgemClusterCol){
217 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
218 return StatusCode::FAILURE;
219 }
220
221 SmartDataPtr<Event::CgemMcHitCol> recCgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
222 if (!recCgemMcHitCol)
223 {
224 log << MSG::WARNING << "Could not find event" << endreq;
225 return StatusCode::FAILURE;
226 }
227
228 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
229 if (!mcParticleCol) {
230 log << MSG::ERROR<< "Could not find McParticle" << endreq;
231 }
232 else{
233
234 Hep3Vector truth_p; HepPoint3D truth_position; double mc_charge(0.);
235 McParticleCol::iterator iter_mc = mcParticleCol->begin();
236 for (;iter_mc != mcParticleCol->end(); ++iter_mc ) {
237 int pdg=(*iter_mc)->particleProperty();
238 if(fabs(pdg)!=13) continue;
239 if(truth_ntrack >=10)continue;
240 if(pdg<0)mc_charge = 1.;
241 else mc_charge =-1.;
242 double truth_x = (*iter_mc)->initialPosition().x();
243 double truth_y = (*iter_mc)->initialPosition().y();
244 double truth_z = (*iter_mc)->initialPosition().z();
245
246 double truth_px = (*iter_mc)->initialFourMomentum().px();
247 double truth_py = (*iter_mc)->initialFourMomentum().py();
248 double truth_pz = (*iter_mc)->initialFourMomentum().pz();
249 double truth_P = sqrt(truth_px*truth_px+truth_py*truth_py+truth_pz*truth_pz);
250 truth_CosTheta[truth_ntrack] = truth_pz/truth_P;
251
252 truth_position = HepPoint3D(truth_x, truth_y, truth_z);
253 truth_p = Hep3Vector(truth_px, truth_py, truth_pz);
254
255 KalmanFit::Helix* tmp_helix = new KalmanFit::Helix(truth_position, truth_p, mc_charge);
256 HepPoint3D tmp_pivot(0., 0., 0.);
257 tmp_helix->pivot(tmp_pivot);
258 truth_dr[truth_ntrack] = tmp_helix->dr();
259 truth_dz[truth_ntrack] = tmp_helix->dz();
260 truth_phi0[truth_ntrack] = tmp_helix->phi0();
261 truth_kappa[truth_ntrack] = tmp_helix->kappa();
262 truth_tanl[truth_ntrack] = tmp_helix->tanl();
263 truth_charge[truth_ntrack] = mc_charge;
264 if(m_debug){
265 cout<<"Track Helix : "<<endl
266 <<"dr = "<<tmp_helix->dr()<<"\tphi0 = "<<tmp_helix->phi0()<<"\tkappa = "<<tmp_helix->kappa()
267 <<"\ndz = "<<tmp_helix->dz()<<"\ttanl = "<<tmp_helix->tanl()<<"\tpt = "<<1./tmp_helix->kappa()<<endl<<endl;
268 }
269
270 truth_ntrack++;
271 delete tmp_helix;
272 }
273 }
274 bool use_truth_hits = true;
275 if(use_truth_hits){
276 SmartDataPtr<Event::CgemMcHitCol> lv_CgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
277 if (!lv_CgemMcHitCol)
278 {
279 log << MSG::WARNING << "Could not find event" << endreq;
280 return StatusCode::FAILURE;
281 }
282 Event::CgemMcHitCol::iterator iter_truth = lv_CgemMcHitCol->begin();
283 for( ; iter_truth != lv_CgemMcHitCol->end(); ++iter_truth){
284 if(fabs((*iter_truth)->GetPDGCode())!=13) continue;
285 double tmp_x = ((*iter_truth)->GetPositionXOfPrePoint()+(*iter_truth)->GetPositionXOfPostPoint())/2;
286 double tmp_y = ((*iter_truth)->GetPositionYOfPrePoint()+(*iter_truth)->GetPositionYOfPostPoint())/2;
287 double tmp_z = ((*iter_truth)->GetPositionZOfPrePoint()+(*iter_truth)->GetPositionZOfPostPoint())/2;
288 double tmp_phi = atan2(tmp_y,tmp_x);
289 while(tmp_phi < 0) tmp_phi += 2*pi;
290 mc_phi[(*iter_truth)->GetLayerID()] = tmp_phi;
291 double tmp_v = ((*iter_truth)->GetPositionZOfPostPoint()+(*iter_truth)->GetPositionZOfPrePoint())/2./10.*sin(CgemSegmentFun::a_stero[(*iter_truth)->GetLayerID()])*CgemSegmentFun::r_layer[(*iter_truth)->GetLayerID()]/CgemSegmentFun::rl[(*iter_truth)->GetLayerID()]+CgemSegmentFun::r_layer[(*iter_truth)->GetLayerID()]*tmp_phi*cos(CgemSegmentFun::a_stero[(*iter_truth)->GetLayerID()]);
292 mc_v[(*iter_truth)->GetLayerID()] = tmp_v;
293 CgemSegmentFun::vec_z[(*iter_truth)->GetLayerID()] = ((*iter_truth)->GetPositionZOfPostPoint()+(*iter_truth)->GetPositionZOfPrePoint())/2./10.;//cm
294 CgemSegmentFun::vec_phi[(*iter_truth)->GetLayerID()] = tmp_phi;
295 // cout<<"truth_z = "<<tmp_z<<"\ttruth_v = "<<tmp_v<<"\ttruth_phi = "<<tmp_phi<<endl;
296 // cout<<"truth_z = "<<tmp_z<<"\ttruth_v = "<<tmp_v<<"\ttruth_phi = "<<tmp_phi<<endl;
297 }
298 double Ini_d0(0.), Ini_phi0(0.), Ini_kappa(0.), Ini_z0(0.), Ini_tanl(0.);
299 CgemSegmentFun::GetHelixVarBeforeFit( -1, Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl);
300 }
301
302 itSeg = recCgemSegmentCol->begin();
303 cgem_nsegment=0;
304 for(;itSeg != recCgemSegmentCol->end(); ++itSeg){
305 if(cgem_nsegment >= 1000){check_segOver++;continue;}
306 cgem_segmentid[cgem_nsegment] = (*itSeg)->getsegmentid();
307 vec_clusterid[0] = (*itSeg)->getclusterid_1();
308 vec_clusterid[1] = (*itSeg)->getclusterid_2();
309 vec_clusterid[2] = (*itSeg)->getclusterid_3();
310 int IsTruth = (*itSeg)->getmatch();
311
312 // initialize minuit
313 Int_t ierflg;
314 Double_t arglist[10];
315 m_failed = 0;
316
317 m_pMnTrk = new TMinuit(5);
318 m_pMnTrk -> SetPrintLevel(-1);
319 m_pMnTrk -> SetFCN(CgemSegmentFun::fcnTrk);
320 m_pMnTrk -> SetErrorDef(1.0); // For Chisq
321 // m_pMnTrk -> SetErrorDef(0.5); // For likelihood
322 m_pMnTrk -> mnparm(0, "dr" , 0, 0.1 , -10, 10, ierflg);//Set starting values and step sizes for parameters
323 m_pMnTrk -> mnparm(1, "phi0" , 0, 0.16, 0, 6.3, ierflg);
324 m_pMnTrk -> mnparm(2, "kappa", 0, 0.4, -100, 100, ierflg);
325 m_pMnTrk -> mnparm(3, "dz" , 0, 0.5 , -50, 50, ierflg);
326 m_pMnTrk -> mnparm(4, "tanL" , 0, 0.005 , -10, 10, ierflg);
327 //m_pMnTrk -> mnparm(0, "dr" , 0, 0.09 , -10, 10, ierflg);//Set starting values and step sizes for parameters
328 //m_pMnTrk -> mnparm(1, "phi0" , 0, 0.016, 0, 6.3, ierflg);
329 //m_pMnTrk -> mnparm(2, "kappa", 0, 0.01, -30, 30, ierflg);
330 //m_pMnTrk -> mnparm(3, "dz" , 0, 0.05 , -50, 50, ierflg);
331 //m_pMnTrk -> mnparm(4, "tanL" , 0, 0.004 , -10, 10, ierflg);
332 arglist[0] = 0;
333 m_pMnTrk -> mnexcm("Set NOW", arglist, 0, ierflg);
334
335
336 RecCgemClusterCol::iterator itTrk=recCgemClusterCol->begin();
337 for(; itTrk != recCgemClusterCol->end(); ++itTrk){
338 if((*itTrk)->getflag() == 0){
339 n_stripX[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
340 CgemSegmentFun::nStrip_X[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;//cm
341 }
342 if((*itTrk)->getflag() == 1){
343 n_stripV[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
344 CgemSegmentFun::nStrip_V[(*itTrk)->getlayerid()] = (*itTrk)->getclusterflage()-(*itTrk)->getclusterflagb()+1;
345 }
346 if((*itTrk)->getflag() != 2) continue;
347 if(!CheckClusterID((*itTrk)->getlayerid(), (*itTrk)->getclusterid(), vec_clusterid)) continue;
348 //CgemSegmentFun::sheet_ID[(*itTrk)->getlayerid()] = (*itTrk)->getsheetid();
349 double tmp_phi = (*itTrk)->getrecphi();
350
351 CgemSegmentFun::vec_z[(*itTrk)->getlayerid()] = (*itTrk)->getRecZ()/10.;//cm
352 CgemSegmentFun::vec_phi[(*itTrk)->getlayerid()] = (*itTrk)->getrecphi();
353 //CgemSegmentFun::vec_v[(*itTrk)->getlayerid()] = (*itTrk)->getrecv()/10.;//cm
354 rec_phi[(*itTrk)->getlayerid()] = (*itTrk)->getrecphi();
355 rec_v[(*itTrk)->getlayerid()] = (*itTrk)->getRecZ()/10.*sin(CgemSegmentFun::a_stero[(*itTrk)->getlayerid()])*CgemSegmentFun::r_layer[(*itTrk)->getlayerid()]/CgemSegmentFun::rl[(*itTrk)->getlayerid()]+CgemSegmentFun::r_layer[(*itTrk)->getlayerid()]*tmp_phi*cos(CgemSegmentFun::a_stero[(*itTrk)->getlayerid()]);//cm
356 CgemSegmentFun::vec_v[(*itTrk)->getlayerid()] = rec_v[(*itTrk)->getlayerid()];//cm
357 //cout<<"vec_Z = "<<(*itTrk)->getRecZ()/10.<<"\tvec_phi = "<<(*itTrk)->getrecphi()<<endl;
358 //cout<<"layerID = "<<(*itTrk)->getlayerid()<<"\tX_n_strip = "<<n_stripX[(*itTrk)->getlayerid()] <<"\tV_n_strip = "<<n_stripV[(*itTrk)->getlayerid()] <<endl;
359 //cout<<"vec_v["<<(*itTrk)->getlayerid()<<"] = "<<(*itTrk)->getrecv()/10.<<"\t"<<rec_v[(*itTrk)->getlayerid()]<<endl;
360 //cout<<"sheetID["<<(*itTrk)->getlayerid()<<"] = "<<(*itTrk)->getsheetid()<<endl;
361 }
362
363 double Ini_d0(0.), Ini_phi0(0.), Ini_kappa(0.), Ini_z0(0.), Ini_tanl(0.);
364 double chisq(0.);
365
366 int charge = GetChargeOfTrack();
367 cgem_charge[cgem_nsegment] = charge;
368 CgemSegmentFun::GetHelixVarBeforeFit( charge, Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl);
369
370 double ErrMatrix[5][5] ;
371 double TrkPar[5] = {Ini_d0, Ini_phi0, Ini_kappa, Ini_z0, Ini_tanl};
372 //double TrkPar[5] = {truth_dr[0], truth_phi0[0], truth_kappa[0], truth_dz[0], truth_tanl[0]};
373 double TrkParErr[5] ;
374 double helixErr[15];
375 int fit = mnTrkFit(TrkPar, TrkParErr, &ErrMatrix[0][0], chisq);
376 int fitAgain = 0; double kappa_err;
377 if(!fit){
378 // cout<< "fit failed"<<"\tevent ID = "<<cgem_evt<<"\trun = "<<cgem_run << endl;
379 m_failed = 1;
380 double tmp_recphi[3], tmp_kappa[8];
381 for(int ii = 0, mm = 0; ii < 2; ii++)
382 for(int jj = 0; jj < 2; jj++)
383 for(int kk = 0; kk < 2; kk++, mm++){
384 tmp_recphi[0] = CgemSegmentFun::vec_phi[0] + pow(-1.,ii)*(wide_delta*n_stripX[0]+pitch*(n_stripX[0]-1))/CgemSegmentFun::r_layer[0];
385 tmp_recphi[1] = CgemSegmentFun::vec_phi[1] + pow(-1.,jj)*(wide_delta*n_stripX[1]+pitch*(n_stripX[1]-1))/CgemSegmentFun::r_layer[1];
386 tmp_recphi[2] = CgemSegmentFun::vec_phi[2] + pow(-1.,kk)*(wide_delta*n_stripX[2]+pitch*(n_stripX[2]-1))/CgemSegmentFun::r_layer[2];
387 charge = GetChargeOfTrack(tmp_recphi);
388 tmp_kappa[mm] = CgemSegmentFun::GetKappaAfterFit(charge, tmp_recphi);
389 //cout<<"tmp kappa["<<mm<<"] = "<<tmp_kappa[mm]<<endl;
390 }
391 double kappa_s = GetMin(tmp_kappa);
392 double kappa_l = GetMax(tmp_kappa);
393 kappa_err = (kappa_s - kappa_l)/2.;
394 double TrkPar[5] = {Ini_d0, Ini_phi0, (kappa_l+kappa_s)/2., Ini_z0, Ini_tanl};
395 TrkPar[2] = kappa_l;
396 fitAgain = mnTrkFit(TrkPar, TrkParErr, &ErrMatrix[0][0], chisq);
397
398 }
399 //cout<<"m_failed = "<<m_failed<<endl;
400
401 //cout<<"charge = "<<charge<<"\t"<<fit<<"\t"<<fitAgain<<endl<<endl;
402 // cout<<
403 // "delta_dr = "<<Ini_d0-truth_dr[0]<<endl<<
404 // "delta_phi0 = "<<Ini_phi0-truth_phi0[0]<<endl<<
405 // "delta_kappa = "<<Ini_kappa-truth_kappa[0]<<endl<<
406 // "delta_dz = "<<Ini_z0-truth_dz[0]<<endl<<
407 // "delta_tanl = "<<Ini_tanl-truth_tanl[0]<<endl<<endl;
408
409 //RecCgemSegment* recsegment = new RecCgemSegment();
410 cgem_fitFlag[cgem_nsegment] = fit;
411 if(fit==1||fitAgain){
412 if(fit){
413 for (int ie = 0, k = 0 ; ie < 5 ; ie++){
414 for (int je = ie ; je < 5 ; je ++, k++) {
415 helixErr[k] = ErrMatrix[ie][je];
416 }
417 }
418 }
419
420 if(fitAgain){
421 for (int ie = 0, k = 0 ; ie < 5 ; ie++){
422 for (int je = ie ; je < 5 ; je ++, k++) {
423 if(ie==2&&je<2)helixErr[k] = ErrMatrix[ie+2][je]; //k = 2,6
424 else if(je<2&&ie>2) helixErr[k] = ErrMatrix[ie-1][je]; //k = 3,4,7,8
425 else if(je>2&&ie>2) helixErr[k] = ErrMatrix[ie-1][je-1]; //k = 12,13,14
426 else if(je==2) helixErr[k] = 0.; //k = 9,10,11
427 else helixErr[k] = ErrMatrix[ie][je]; // k = 0,1,5
428
429 if(k==9)helixErr[k] = pow(kappa_err,2);
430 // cout<<"ErrMatrix["<<ie<<"]["<<je<<"] = "<<ErrMatrix[ie][je]<<"\t"<<helixErr[k]<<endl;
431 }
432 }
433 TrkParErr[0] = 0.8812/0.0677*TrkParErr[0];
434 TrkParErr[1] = 0.266/0.004*TrkParErr[1];
435 TrkParErr[2] = kappa_err;
436 helixErr[0] = 0.8812/0.0677*helixErr[0];
437 helixErr[5] = 0.266/0.004*helixErr[5];
438 }
439 // cout<<"helixErr[0] = "<<helixErr[0]<<"\tTrkParErr[0] = "<<TrkParErr[0]<<"\n"<<
440 // "helixErr[5] = "<<helixErr[5]<<"\tTrkParErr[1] = "<<TrkParErr[1]<<"\n"<<
441 // "helixErr[9] = "<<helixErr[9]<<"\tTrkParErr[2] = "<<TrkParErr[2]<<"\n"<<
442 // "helixErr[12] = "<<helixErr[12]<<"\tTrkParErr[3] = "<<TrkParErr[3]<<"\n"<<
443 // "helixErr[14] = "<<helixErr[14]<<"\tTrkParErr[4] = "<<TrkParErr[4]<<"\n"<<endl<<endl;
444 //recsegment->sethelix(TrkPar);
445 //recsegment->sethelix_err(helixErr);
446 //recsegment->setmatch(IsTruth);
447 (*itSeg)->sethelix(TrkPar);
448 (*itSeg)->sethelix_err(helixErr);
449 (*itSeg)->setmatch(IsTruth);
450 cgem_dr[cgem_nsegment] = TrkPar[0];
451 cgem_phi0[cgem_nsegment] = TrkPar[1];
452 cgem_kappa[cgem_nsegment] = TrkPar[2];
453 cgem_dz[cgem_nsegment] = TrkPar[3];
454 cgem_tanl[cgem_nsegment] = TrkPar[4];
455 cgem_chisq[cgem_nsegment] = chisq;
456 if(m_debug){
457 cout<< "fit success" << endl
458 <<"dr = " <<TrkPar[0]<<" +- "<<TrkParErr[0]<<"\t\t"
459 <<"phi0 = "<<TrkPar[1]<<" +- "<<TrkParErr[1]<<"\t\t"
460 <<"kappa= "<<TrkPar[2]<<" +- "<<TrkParErr[2]<<"\n"
461 <<"dz = " <<TrkPar[3]<<" +- "<<TrkParErr[3]<<"\t\t"
462 <<"tanl = "<<TrkPar[4]<<" +- "<<TrkParErr[4]<<"\n\n"
463 <<"kappa_err= "<<TrkParErr[2]<<endl;
464 cout<<"chisq = "<<chisq<<endl;
465 }
466 }
467 else{
468 //cout<< "fit failed"<<"\tevent ID = "<<cgem_evt<<"\trun = "<<cgem_run << endl;
469
470 double fail_helix[5], fail_helix_err[15] ;
471 for(int i = 0; i <15; i++){
472 if(i<5)fail_helix[i] = -999.;
473 fail_helix_err[i] = -999.;
474 }
475 (*itSeg)->sethelix(fail_helix);
476 (*itSeg)->sethelix_err(fail_helix_err);
477 (*itSeg)->setmatch(IsTruth);
478 cgem_dr[cgem_nsegment] = -999.;
479 cgem_phi0[cgem_nsegment] = -999.;
480 cgem_kappa[cgem_nsegment] = -999.;
481 cgem_dz[cgem_nsegment] = -999.;
482 cgem_tanl[cgem_nsegment] = -999.;
483 cgem_chisq[cgem_nsegment] = -999.;
484 }
485 //m_seg_map.insert(pair<int,RecCgemSegment*>(cgem_nsegment,recsegment));
486
487 cgem_nsegment++;
488
489 //cout<<"***************************#######*******************************"<<endl<<endl;
490 delete m_pMnTrk;
491 //delete recsegment;
492 }
493 /*
494 RecCgemSegmentCol* segmentcol = new RecCgemSegmentCol;
495
496 IDataProviderSvc* evtSvc = NULL;
497 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
498 if (evtSvc) {
499 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
500 } else {
501 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
502 return StatusCode::SUCCESS;
503 }
504
505 StatusCode segmentsc;
506 IDataManagerSvc *dataManSvc;
507 dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
508 DataObject *aSegmentCol;
509 evtSvc->findObject("/Event/Recon/RecCgemSegmentCol",aSegmentCol);
510 if(aSegmentCol != NULL) {
511 dataManSvc->clearSubTree("/Event/Recon/RecCgemSegmentCol");
512 evtSvc->unregisterObject("/Event/Recon/RecCgemSegmentCol");
513 }
514
515
516 //push back the segment from the map//
517 for(m_seg_map_it = m_seg_map.begin();m_seg_map_it!=m_seg_map.end();++m_seg_map_it){
518 RecCgemSegment* recsegment = (*m_seg_map_it).second;
519 segmentcol->push_back(recsegment);
520 cout<<">>>>>>> push recsegment, clusterid_1="<<recsegment->getclusterid_1()<<"kappa: "<<recsegment->gethelix(2)<<endl;
521
522 }
523
524 segmentsc = evtSvc->registerObject("/Event/Recon/RecCgemSegmentCol", segmentcol);
525 if( segmentsc.isFailure() ) {
526 log << MSG::FATAL << "Could not register RecCgemSegment" << endreq;
527 return StatusCode::SUCCESS;
528 }
529 log << MSG::INFO << "RecCgemSegmentCol registered successfully!" <<endreq;
530 */
531
532
533
534 helix_ntuple->write();
535 /*
536 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol2(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
537 itSeg = recCgemSegmentCol2->begin();
538 for(;itSeg != recCgemSegmentCol2->end(); itSeg++){
539 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
540 cout<<"event ID = "<<(*itSeg)->geteventid()<<endl;;
541 cout<<"segment ID = "<<(*itSeg)->getsegmentid()<<endl;;
542 cout<<"cluster ID_1 = "<<(*itSeg)->getclusterid_1()<<endl;;
543 cout<<"cluster ID_2 = "<<(*itSeg)->getclusterid_2()<<endl;;
544 cout<<"cluster ID_3 = "<<(*itSeg)->getclusterid_3()<<endl;;
545 cout<<"event+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
546 }
547 */
548
549
550 return StatusCode::SUCCESS;
551}
552
553StatusCode CgemSegmentFitAlg::exe_v2()
554{
555 MsgStream log(msgSvc(), name());
556 log << MSG::INFO << "in excute()" << endreq;
557
558 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
559 if (!eventHeader) {
560 log << MSG::FATAL << "Could not find Event Header" << endreq;
561 return StatusCode::FAILURE;
562 }
563 m_run = eventHeader->runNumber();
564 m_evt = eventHeader->eventNumber();
565 getMcPar();
566
567 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
568 if (!recCgemClusterCol){
569 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
570 return StatusCode::FAILURE;
571 }
572 RecCgemClusterCol::iterator itCluBegin=recCgemClusterCol->begin();
573
574 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
575 if (!recCgemSegmentCol){
576 log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
577 check_segNum++;
578 // return StatusCode::FAILURE;
579 return StatusCode::SUCCESS;
580 }
581 if(m_checkIdx) {
582 cgem_run = m_run;
583 cgem_evt = m_evt;
584 cgem_nsegment=0;
585 }
586 int cluId[3];
587 RecCgemSegmentCol::iterator itSeg = recCgemSegmentCol->begin();
588 for( ;itSeg != recCgemSegmentCol->end(); ++itSeg)
589 {
590 cluId[0] = (*itSeg)->getclusterid_1();
591 cluId[1] = (*itSeg)->getclusterid_2();
592 cluId[2] = (*itSeg)->getclusterid_3();
593 if(m_debug)
594 {
595 cout<<"input clusters: "<<endl;
596 }
597 double vec_Q[3], vec_T[3];
598 for(int i=0; i<3; i++)
599 {
600 RecCgemClusterCol::iterator itClu = itCluBegin+cluId[i];
601 CgemSegmentFun::vec_phi[i] = (*itClu)->getrecphi();// in radian
602 CgemSegmentFun::vec_z[i] = (*itClu)->getRecZ()/10.; // in cm
603 vec_Q[i]=(*itClu)->getenergydeposit();
604 vec_T[i]=100;
605 if(m_debug)
606 {
607 cout <<"layer "<<i<<", cluster id = "<<cluId[i]<<", itCluId = "<<(*itClu)->getclusterid()
608 <<", phi="<<CgemSegmentFun::vec_phi[i]<<", z="<<CgemSegmentFun::vec_z[i]<<endl;
609 }
610 }
611 //cout<<"charge = "<<GetChargeOfTrack()<<", "<<estiChargeOfSeg()<<endl;
612 int charge = estiChargeOfSeg();
613
614 // initial Helix parameters
615 double dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini;
616 CgemSegmentFun::GetHelixVarBeforeFit(charge, dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini);
617 double dr_ini0 = dr_ini;
618 double phi0_ini0 = phi0_ini;
619 double kappa_ini0 = kappa_ini;
620 double dZ_ini0 = dZ_ini;
621 double tanL_ini0 = tanL_ini;
622
623 double helix[5] = {dr_ini, phi0_ini, kappa_ini, dZ_ini, tanL_ini};
624 double helixErr[15];
625 HepPoint3D pivot(0., 0., 0.);
626 HepVector vec_trkpar(5);
627 for(int i = 0; i < 5; i++){
628 vec_trkpar[i] = helix[i];
629 }
630 //KalmanFit::Helix* cgem_helix = new KalmanFit::Helix(pivot, vec_trkpar);
631 KalmanFit::Helix cgem_helix(pivot, vec_trkpar);
632
633 /* initialize error matrix */
634 int iView(0), mode(2);
635 for(int i=0; i<3; i++)
636 {
637 double dphi = CgemSegmentFun::IntersectCylinder(CgemSegmentFun::rl[i], cgem_helix);
638 HepPoint3D pos = cgem_helix.x(dphi);
639 double phi_pos = pos.phi();
640 Hep3Vector p3 = cgem_helix.momentum(dphi);
641 double phi_p = p3.phi();
642 double dPhi = phi_p-phi_pos;
643 while(dPhi>CLHEP::pi) dPhi-=CLHEP::twopi;
644 while(dPhi<-CLHEP::pi) dPhi+=CLHEP::twopi;
645 CgemSegmentFun::x_reso[i]= (myCgemCalibSvc->getSigma(i,0,mode,dPhi,vec_Q[i],vec_T[i]))/10.;// mm -> cm
646 CgemSegmentFun::v_reso[i]= (myCgemCalibSvc->getSigma(i,1,mode,dPhi,vec_Q[i],vec_T[i]))/10.;// mm -> cm
647 }
648
649 // initialize minuit
650 Int_t ierflg, istat;
651 Double_t arglist[10];
652 m_failed = 0;
653
654 m_pMnTrk = new TMinuit(5);// num of par
655 m_pMnTrk -> SetPrintLevel(-1);// no print out
656 m_pMnTrk -> mnexcm("Set NOW", arglist, 1, ierflg);// no warning
657 m_pMnTrk -> SetFCN(CgemSegmentFun::fcnTrk);
658 m_pMnTrk -> SetErrorDef(1.0); // For Chisq
659 // m_pMnTrk -> SetErrorDef(0.5); // For likelihood
660 bool fitSucceed = false;
661 bool fitUsable = false;
662 int istat_last = 0;
663 int nFitTried = 0; int iFitStored = 0;
664 double chi2_min = 999;
665 double chi2_saved = 999;
666 double chi2_fixKapa_min = 990;
667 while((nFitTried<3&&!fitSucceed)||(nFitTried<5&&!fitUsable)) {
668 m_pMnTrk-> mnparm(0, "dr" , dr_ini, 0.001 , dr_ini-10, dr_ini+10, ierflg);//Set starting values and step sizes for parameters
669 m_pMnTrk-> mnparm(1, "phi0" , phi0_ini, 0.001, phi0_ini-1.0, phi0_ini+1.0, ierflg);
670 m_pMnTrk-> mnparm(2, "kappa", kappa_ini, max(fabs(kappa_ini)*0.01,0.01), kappa_ini-5.0, kappa_ini+5.0,ierflg);
671 if(nFitTried>=3) m_pMnTrk->FixParameter(2);
672 m_pMnTrk -> mnparm(3, "dz" , dZ_ini, 0.01 , dZ_ini-20.0, dZ_ini+20.0, ierflg);
673 m_pMnTrk -> mnparm(4, "tanL" , tanL_ini, 0.001 , tanL_ini-1.0, tanL_ini+1.0, ierflg);
674 arglist[0] = 2000; //arglist[1] = 0.1;
675 m_pMnTrk -> mnexcm("MIGRAD", arglist, 1, ierflg);
676 nFitTried++;
677 double trkpar[5], trkparErr[5], emat[25], chisq;
678 Double_t fmin, edm, errdef;
679 Int_t nvpar, nparx;
680 m_pMnTrk -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
681 for(int i=0; i<5; i++){
682 m_pMnTrk -> GetParameter(i, trkpar[i], trkparErr[i]);
683 }
684 // error matrix
685 m_pMnTrk -> mnemat(emat, 5);
686 chisq = fmin;
687
688 //HepVector aNew = CgemSegmentFun::CalculateHelix(vec_trkpar, pivot);
689
690 if(m_debug) {
691 cout<<"Fit ierflg="<<ierflg<<", istat="<<istat<<endl;
692 cout.precision(4);
693 cout<<setw(10)<<" "<<setw(10)<<"dr"<<setw(10)<<"phi0"<<setw(10)<<"kappa"<<setw(10)<<"dz"<<setw(10)<<"tanL"<<endl;
694 cout<<setw(10)<<"MC truth"<<setw(10)<<m_helix_mc[0]<<setw(10)<<m_helix_mc[1]<<setw(10)<<m_helix_mc[2]<<setw(10)<<m_helix_mc[3]<<setw(10)<<m_helix_mc[4]<<endl;
695 cout<<setw(10)<<"ini"<<setw(10)<<dr_ini<<setw(10)<<phi0_ini<<setw(10)<<kappa_ini<<setw(10)<<dZ_ini<<setw(10)<<tanL_ini<<endl;
696 //cout<<setw(10)<<"calcu"<<setw(10)<<aNew[0]<<setw(10)<<aNew[1]<<setw(10)<<aNew[2]<<setw(10)<<aNew[3]<<setw(10)<<aNew[4]<<"chi2="<<chisq<<endl;
697 cout<<setw(10)<<"fit"<<setw(10)<<trkpar[0]<<setw(10)<<trkpar[1]<<setw(10)<<trkpar[2]<<setw(10)<<trkpar[3]<<setw(10)<<trkpar[4]<<"chi2="<<chisq<<endl;
698 cout<<setw(10)<<"fit err"<<setw(10)<<trkparErr[0]<<setw(10)<<trkparErr[1]<<setw(10)<<trkparErr[2]<<setw(10)<<trkparErr[3]<<setw(10)<<trkparErr[4]<<endl;
699 cout.precision(6);
700 }
701
702 /*
703 *-* IERFLG is now (94.5) defined the same as ICONDN in MNCOMD
704 *-* = 0: command executed normally
705 *-* 1: command is blank, ignored
706 *-* 2: command line unreadable, ignored
707 *-* 3: unknown command, ignored
708 *-* 4: abnormal termination (e.g., MIGRAD not converged)
709 *-* 9: reserved
710 *-* 10: END command
711 *-* 11: EXIT or STOP command
712 *-* 12: RETURN command
713 *
714 * ISTAT: a status integer indicating how good is the covariance
715 *-* matrix: 0= not calculated at all
716 *-* 1= approximation only, not accurate
717 *-* 2= full matrix, but forced positive-definite
718 *-* 3= full accurate covariance matrix
719 */
720 if( (0==ierflg) && (istat==3) ){
721 fitSucceed=true;
722 }
723 else {
724 //cout<<"TMinuit failed"<<endl;
725 if(m_debug) {
726 cout<<"TMinuit tried "<<nFitTried<<" times!"<<endl;
727 }
728 dr_ini=trkpar[0];
729 phi0_ini=trkpar[1];
730 //kappa_ini=trkpar[2]+pow(-1,nFitTried);
731 kappa_ini=trkpar[2]+0.5*pow(-1,nFitTried);
732 dZ_ini=trkpar[3];
733 tanL_ini=trkpar[4];
734 }
735 if(nFitTried>=3) kappa_ini=0.5*pow(-1,nFitTried);
736
737 if( (0==ierflg) && (istat>0) ){
738 if(nFitTried<4) fitUsable = true;
739 if((nFitTried<4&&istat>istat_last)||(nFitTried>3&&chisq<chi2_fixKapa_min)) {
740 int k=0;
741 chi2_fixKapa_min = chisq;
742 chi2_saved = chisq;
743 for(int ii=0; ii<5; ii++)
744 {
745 helix[ii]=trkpar[ii];
746 for(int jj=ii; jj<5; jj++)
747 {
748 helixErr[k]=emat[ii*5+jj];
749 k++;
750 }
751 }
752 istat_last = istat;
753 iFitStored = nFitTried;
754
755 }
756 }
757 }
758 //if(!fitUsable) cout<<"TMinuit failed after "<<nFitTried<<" attemps!"<<endl;
759 if(m_debug)
760 cout<<"Fit results "<<iFitStored<<" are saved!"<<endl;
761
762 if(m_checkIdx&&cgem_nsegment<1000) {
763 cgem_segmentid[cgem_nsegment]=(*itSeg)->getsegmentid();
764 cgem_dr_ini[cgem_nsegment]=dr_ini0;
765 cgem_dz_ini[cgem_nsegment]=dZ_ini0;
766 cgem_phi0_ini[cgem_nsegment]=phi0_ini0;
767 cgem_kappa_ini[cgem_nsegment]=kappa_ini0;
768 cgem_tanl_ini[cgem_nsegment]=tanL_ini0;
769 cgem_dr[cgem_nsegment]=helix[0];
770 cgem_dz[cgem_nsegment]=helix[3];
771 cgem_phi0[cgem_nsegment]=helix[1];
772 cgem_kappa[cgem_nsegment]=helix[2];
773 cgem_tanl[cgem_nsegment]=helix[4];
774 cgem_chisq[cgem_nsegment]=chi2_saved;
775 cgem_fitFlag[cgem_nsegment]=iFitStored;
776 cgem_charge[cgem_nsegment]=charge;
777 cgem_nsegment++;
778 }
779
780 // fill out
781 (*itSeg)->sethelix(helix);
782 (*itSeg)->sethelix_err(helixErr);
783
784 delete m_pMnTrk;
785
786 }// loop CGEM segments
787 //cout<<"end exe_v2()"<<endl;
788 if(m_checkIdx) helix_ntuple->write();
789
790 return StatusCode::SUCCESS;
791
792}
793
794double CgemSegmentFitAlg::GetMin(double *a){
795 double tmp = 999.;
796 for(int i = 0; i < 8; i++){
797 if(a[i] <= tmp) tmp = a[i];
798 }
799 return tmp;
800}
801
802double CgemSegmentFitAlg::GetMax(double *a){
803 double tmp = -999.;
804 for(int i = 0; i < 8; i++){
805 if(a[i] >= tmp) tmp = a[i];
806 }
807 return tmp;
808}
809
810bool CgemSegmentFitAlg::CheckClusterID(int layerID, int clusterID, int vec_clusterID[3]){
811 bool IsEq = false;
812 if(clusterID == vec_clusterID[layerID]) IsEq = true;
813 return IsEq;
814}
815
816int CgemSegmentFitAlg::GetChargeOfTrack(double *vec_phi){
817 bool IsCross = false;
818 double phi31 = atan2(CgemSegmentFun::rl[2]*sin(vec_phi[2])-CgemSegmentFun::rl[0]*sin(vec_phi[0]), CgemSegmentFun::rl[2]*cos(vec_phi[2])-CgemSegmentFun::rl[0]*cos(vec_phi[0]));
819 double phi21 = atan2(CgemSegmentFun::rl[1]*sin(vec_phi[1])-CgemSegmentFun::rl[0]*sin(vec_phi[0]), CgemSegmentFun::rl[1]*cos(vec_phi[1])-CgemSegmentFun::rl[0]*cos(vec_phi[0]));
820 double phi32 = atan2(CgemSegmentFun::rl[2]*sin(vec_phi[2])-CgemSegmentFun::rl[1]*sin(vec_phi[1]), CgemSegmentFun::rl[2]*cos(vec_phi[2])-CgemSegmentFun::rl[1]*cos(vec_phi[1]));
821 if(phi31*phi21 < 0.) IsCross = true;
822 if(phi31 < 0.)phi31 += 2.*M_PI;
823 if(phi21 < 0.)phi21 += 2.*M_PI;
824 if(phi32 < 0.)phi32 += 2.*M_PI;
825 if(IsCross && phi21 >= phi32 && phi32 >= phi31) return 1;
826 if(IsCross && phi31 >= phi32 && phi32 >= phi21) return -1;
827 if(IsCross && phi21 >= phi31 && phi31 >= phi32) return -1;
828 if(IsCross && phi32 >= phi31 && phi31 >= phi21) return 1;
829 if(!(IsCross) && phi31 >= phi21) return 1;
830 if(!(IsCross) && phi21 > phi31) return -1;
831}
832int CgemSegmentFitAlg::GetChargeOfTrack(){
833 bool IsCross = false;
834 double phi31 = atan2(CgemSegmentFun::rl[2]*sin(CgemSegmentFun::vec_phi[2])-CgemSegmentFun::rl[0]*sin(CgemSegmentFun::vec_phi[0]), CgemSegmentFun::rl[2]*cos(CgemSegmentFun::vec_phi[2])-CgemSegmentFun::rl[0]*cos(CgemSegmentFun::vec_phi[0]));
835 double phi21 = atan2(CgemSegmentFun::rl[1]*sin(CgemSegmentFun::vec_phi[1])-CgemSegmentFun::rl[0]*sin(CgemSegmentFun::vec_phi[0]), CgemSegmentFun::rl[1]*cos(CgemSegmentFun::vec_phi[1])-CgemSegmentFun::rl[0]*cos(CgemSegmentFun::vec_phi[0]));
836 double phi32 = atan2(CgemSegmentFun::rl[2]*sin(CgemSegmentFun::vec_phi[2])-CgemSegmentFun::rl[1]*sin(CgemSegmentFun::vec_phi[1]), CgemSegmentFun::rl[2]*cos(CgemSegmentFun::vec_phi[2])-CgemSegmentFun::rl[1]*cos(CgemSegmentFun::vec_phi[1]));
837 if(phi31*phi21 < 0.) IsCross = true;
838 if(phi31 < 0.)phi31 += 2.*M_PI;
839 if(phi21 < 0.)phi21 += 2.*M_PI;
840 if(phi32 < 0.)phi32 += 2.*M_PI;
841 if(IsCross && phi21 >= phi32 && phi32 >= phi31) return 1;
842 if(IsCross && phi31 >= phi32 && phi32 >= phi21) return -1;
843 if(IsCross && phi21 >= phi31 && phi31 >= phi32) return -1;
844 if(IsCross && phi32 >= phi31 && phi31 >= phi21) return 1;
845 if(!(IsCross) && phi31 >= phi21) return 1;
846 if(!(IsCross) && phi21 > phi31) return -1;
847}
848
849int CgemSegmentFitAlg::estiChargeOfSeg()
850{
851 double x[3],y[3];
852 for(int i=0; i<3; i++)
853 {
854 x[i]=CgemSegmentFun::rl[i]*cos(CgemSegmentFun::vec_phi[i]);
855 y[i]=CgemSegmentFun::rl[i]*sin(CgemSegmentFun::vec_phi[i]);
856 }
857 double x21=x[1]-x[0];
858 double y21=y[1]-y[0];
859 double phi21=atan2(y21, x21);
860 double x32=x[2]-x[1];
861 double y32=y[2]-y[1];
862 double phi32=atan2(y32, x32);
863 double dPhi = phi32-phi21;
864 while(dPhi>acos(-1.)) dPhi-=2.0*acos(-1.);
865 while(dPhi<-acos(-1.)) dPhi+=2.0*acos(-1.);
866 if(dPhi==0.0) return 0;
867 else return dPhi/fabs(dPhi);
868}
869
870
871int CgemSegmentFitAlg::mnTrkFit(Double_t trkpar[], Double_t trkparErr[],
872 Double_t* emat, double& chisq){
873
874 Int_t ierflg;
875 Int_t istat;
876 Double_t arglist[10];
877 Double_t fmin;
878 Double_t edm;
879 Double_t errdef;
880 Int_t nvpar;
881 Int_t nparx;
882 // int i =0;
883 // while(i<5)cout<<"init trk par["<<i<<"] ="<<trkpar[i++]<<endl;
884
885 arglist[0] = 1;
886 arglist[1] = trkpar[0];
887 m_pMnTrk -> mnexcm("SET PARameter", arglist, 2, ierflg);
888 arglist[0] = 2;
889 arglist[1] = trkpar[1];
890 m_pMnTrk -> mnexcm("SET PARameter", arglist, 2, ierflg);
891 arglist[0] = 3;
892 arglist[1] = trkpar[2];
893 m_pMnTrk -> mnexcm("SET PARameter", arglist, 2, ierflg);
894 if(m_failed) m_pMnTrk->FixParameter(2);
895 //if(m_failed) m_pMnTrk->mnexcm("FIX PARameter", arglist, 2, ierflg);
896 arglist[0] = 4;
897 arglist[1] = trkpar[3];
898 m_pMnTrk -> mnexcm("SET PARameter", arglist, 2, ierflg);
899 arglist[0] = 5;
900 arglist[1] = trkpar[4];
901 m_pMnTrk -> mnexcm("SET PARameter", arglist, 2, ierflg);
902
903 //if(m_failed) m_pMnTrk->FixParameter(3);
904 arglist[0] = 2000;
905 arglist[1] = 0.1;
906 m_pMnTrk -> mnexcm("MIGRAD", arglist, 2, ierflg);
907 // if(m_failed){
908 // m_pMnTrk->Release(2);
909 // // m_pMnTrk -> mnscan();
910 // arglist[0] = 3;
911 // arglist[1] = trkpar[2];
912 // m_pMnTrk -> mnexcm("SCAN", arglist, 1, ierflg);
913 // }
914 //m_pMnTrk -> mnexcm("MINIMIZE", arglist, 2, ierflg);
915 m_pMnTrk -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
916
917 //cout<<"ierrflg = "<<ierflg<<"\tistat = "<<istat<<endl;
918
919 // if( (0==ierflg) && (2==istat) ){
920 // for(int i=0; i<5; i++){
921 // m_pMnTrk -> GetParameter(i, trkpar[i], trkparErr[i]);
922 // }
923 //
924 // // error matrix
925 // m_pMnTrk -> mnemat(emat, 5);
926 // chisq = fmin;
927 //
928 // return 2;
929 // }
930 if( (0==ierflg) && (3==istat) ){
931 for(int i=0; i<5; i++){
932 m_pMnTrk -> GetParameter(i, trkpar[i], trkparErr[i]);
933 }
934
935 // error matrix
936 m_pMnTrk -> mnemat(emat, 5);
937 chisq = fmin;
938
939 return 1;
940 }
941
942 return 0;
943}
944
945void CgemSegmentFitAlg::getMcPar()
946{
947 m_helix_mc[2]=99.0;
948 if(m_checkIdx) truth_ntrack = 0;
949 if(m_run<0)
950 {
951 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
952 if (mcParticleCol) {
953 Hep3Vector truth_p; HepPoint3D truth_position; double mc_charge(0.);
954 McParticleCol::iterator iter_mc = mcParticleCol->begin();
955 for (;iter_mc != mcParticleCol->end(); ++iter_mc ) {
956 int pdg=(*iter_mc)->particleProperty();
957 if(fabs(pdg)!=13) continue;
958 if(pdg<0)mc_charge = 1.;
959 else mc_charge =-1.;
960 double truth_x = (*iter_mc)->initialPosition().x();// cm
961 double truth_y = (*iter_mc)->initialPosition().y();
962 double truth_z = (*iter_mc)->initialPosition().z();
963
964 double truth_px = (*iter_mc)->initialFourMomentum().px();
965 double truth_py = (*iter_mc)->initialFourMomentum().py();
966 double truth_pz = (*iter_mc)->initialFourMomentum().pz();
967 double truth_P = sqrt(truth_px*truth_px+truth_py*truth_py+truth_pz*truth_pz);
968 truth_position = HepPoint3D(truth_x, truth_y, truth_z);
969 truth_p = Hep3Vector(truth_px, truth_py, truth_pz);
970 KalmanFit::Helix* tmp_helix = new KalmanFit::Helix(truth_position, truth_p, mc_charge);
971 HepPoint3D tmp_pivot(0., 0., 0.);
972 tmp_helix->pivot(tmp_pivot);
973 m_helix_mc[0] = tmp_helix->dr();
974 m_helix_mc[1] = tmp_helix->phi0();
975 m_helix_mc[2] = tmp_helix->kappa();
976 m_helix_mc[3] = tmp_helix->dz();
977 m_helix_mc[4] = tmp_helix->tanl();
978 if(m_checkIdx&&truth_ntrack<10) {
979 truth_dr[truth_ntrack]=m_helix_mc[0];
980 truth_dz[truth_ntrack] = m_helix_mc[3];
981 truth_phi0[truth_ntrack] = m_helix_mc[1];
982 truth_kappa[truth_ntrack] = m_helix_mc[2];
983 truth_tanl[truth_ntrack] = m_helix_mc[4];
984 truth_charge[truth_ntrack] = mc_charge;
985 truth_CosTheta[truth_ntrack] = truth_pz/truth_P;
986 truth_ntrack++;
987 }
988 if(m_debug){
989 cout<<"Track Helix : "<<endl
990 <<"dr = "<<tmp_helix->dr()<<"\tphi0 = "<<tmp_helix->phi0()<<"\tkappa = "<<tmp_helix->kappa()
991 <<"\ndz = "<<tmp_helix->dz()<<"\ttanl = "<<tmp_helix->tanl()<<"\tpt = "<<1./tmp_helix->kappa()<<endl<<endl;
992 }
993 delete tmp_helix;
994 }
995 }
996 }
997}
998
999// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1000StatusCode CgemSegmentFitAlg::finalize() {
1001
1002 MsgStream log(msgSvc(), name());
1003 log << MSG::INFO << "in finalize()" << endreq;
1004 //cout<<"leave segment number by dZ dPhi = "<<check_segNum<<endl
1005 // <<"leave segment number over 1000 = "<<check_segOver<<endl;
1006
1007 //cout<<endl<<"******* CgemSegmentFitAlg *******"<<endl;
1008
1009 return StatusCode::SUCCESS;
1010}
1011
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Gam4pikp.cxx:37
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
M_PI
#define M_PI
Definition TConstant.h:4
CgemGeoLayer::getMiddleROfGapD
double getMiddleROfGapD() const
Definition CgemGeoLayer.h:105
CgemGeoLayer::getAngleOfStereo
double getAngleOfStereo() const
Definition CgemGeoLayer.h:218
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition CgemGeomSvc.h:48
CgemGeomSvc::getNumberOfCgemLayer
int getNumberOfCgemLayer() const
Definition CgemGeomSvc.h:45
CgemGeomSvc::getReadoutPlane
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition CgemGeomSvc.h:140
CgemSegmentFitAlg::CgemSegmentFitAlg
CgemSegmentFitAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition CgemSegmentFitAlg.cxx:64
ICgemCalibFunSvc::getSigma
virtual double getSigma(int layer, int xvFlag, int readoutMode, double angle, double Q, double T) const =0
KalmanFit::Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-18/src/helix/Helix.cxx:206
KalmanFit::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-18/src/helix/Helix.cxx:151
KalmanFit::Helix::dr
double dr(void) const
returns an element of parameters.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-18/KalFitAlg/helix/Helix.h:229
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-18/KalFitAlg/helix/Helix.h:211
CgemSegmentFun::GetHelixVarBeforeFit
void GetHelixVarBeforeFit(int charge, double &d0, double &phi0, double &omega, double &z0, double &tanl)
Definition CgemSegmentFun.cxx:255
CgemSegmentFun::GetKappaAfterFit
double GetKappaAfterFit(int charge, double *rec_phi)
Definition CgemSegmentFun.cxx:240
CgemSegmentFun::fcnTrk
void fcnTrk(int &npar, double *gin, double &f, double *par, int iflag)
Definition CgemSegmentFun.cxx:205
CgemSegmentFun::IntersectCylinder
double IntersectCylinder(double r, KalmanFit::Helix helix)
Definition CgemSegmentFun.cxx:225
Event
Definition Event.h:21
pi
const float pi
Definition vector3.h:133