CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemSegmentRecAlg.cxx
Go to the documentation of this file.
1#include "CgemSegmentRecAlg/CgemSegmentRecAlg.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 "CLHEP/Vector/ThreeVector.h"
27
28#include "CgemRecEvent/RecCgemCluster.h"
29#include "CgemRecEvent/RecCgemSegment.h"
30//#include "CgemClusterCreat/CgemClusterCreat.h"
31#include "McTruth/CgemMcHit.h"
32#include "CgemRawEvent/CgemDigi.h"
33#include "CgemGeomSvc/CgemGeomSvc.h"
34
35
36#include "VertexFit/KalmanKinematicFit.h"
37#include "VertexFit/VertexFit.h"
38#include "VertexFit/Helix.h"
39
40#include "TRandom.h"
41#include "TArrayI.h"
42
43#include <vector>
44#include <iostream>
45#include <fstream>
46#include <iomanip>
47#include <cstring>
48
49using namespace std;
50using namespace Event;
51
52//vector<int> continuousHit[43];
53
54/////////////////////////////////////////////////////////////////////////////
55CgemSegmentRecAlg::CgemSegmentRecAlg(const std::string& name, ISvcLocator* pSvcLocator) :
56 Algorithm(name, pSvcLocator)
57{
58 declareProperty("Debug", m_debug =0);
59 declareProperty("version", m_ver =2);
60 declareProperty("CheckMCSoleTrk", m_checkSingleTrkMC=0);
61 declareProperty("HistFile", m_fname);
62 declareProperty("HistFileMode", m_fmode = "recreate");
63 declareProperty("PostFix", m_postfix = "");
64 declareProperty("ParticleType", m_pid);
65 declareProperty("isDst", m_dst = true );
66 declareProperty("UseCut", m_UseCut );
67 declareProperty("z_slope", m_z_slope );
68 declareProperty("z_Par1", m_z_Par1 );
69 declareProperty("z_Par2", m_z_Par2 );
70 declareProperty("z_Par3", m_z_Par3 );
71 declareProperty("z_nSigma",m_z_nSigma );
72 declareProperty("z_bound", m_z_bound );
73 declareProperty("phi_slope", m_phi_slope );
74 declareProperty("phi_3Par", m_phi_3Par );
75 declareProperty("phi_4Par", m_phi_4Par );
76 declareProperty("phi_nSigma", m_phi_nSigma );
77 declareProperty("phi_bound", m_phi_bound );
78 declareProperty("UseTesterSim", m_UseTesterSim );
79}
80
81// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
82StatusCode CgemSegmentRecAlg::initialize(){
83
84 MsgStream log(msgSvc(), name());
85 log << MSG::INFO << "in initialize()" << endreq;
86 StatusCode sc;
87
88 Ncut0 = 0;
89 cgem_totevt = 0;
90 N_cluster = 0;N_strip = 0;N_XV_cluster = 0; N_hit = 0;N_flag3 = 0; N_sim_strip = 0; N_trk = 0;N_N21 = 0; N_N23 = 0; N_sim_N21 =0; N_sim_N23 =0;N_OutRange = 0;
91 N_PiL2 = 0;N_PiE0 = 0;N_PiM7 = 0;N_Match=0;N_Event=0;N_4cluster = 0;
92 N_nonPi = 0; N_ErrMat = 0; N_CrossMat = 0; N_OtherPar4 = 0;N_OtherPar3 = 0; N_OtherPar = 0; N_3cluster = 0; N_inTrk3 = 0; N_InTrk = 0; N_outTrk = 0;
93 N_Match_hit = 0;N_Match_all_hit = 0;OutDetecter = 0;
94 if(m_checkSingleTrkMC>0)
95 {
96 /////////////////////////////// book ntuple //////////////////////////////////////
97 INTupleSvc* ntupleSvc;
98 sc = service("NTupleSvc", ntupleSvc);
99 int postfixlen = m_postfix.length();
100 // free(foldername);
101 char *foldername = (char *)malloc((postfixlen + 15) * sizeof(char));
102 sprintf(foldername, "FILE145/CgemRec%s", m_postfix.c_str());
103 NTuplePtr mt(ntupleSvc, foldername);
104 if ( mt ) cgem_ntuple = mt;
105 else{
106 cgem_ntuple = ntupleSvc->book(foldername, CLID_ColumnWiseTuple, "CgemSegment");
107 // NTuplePtr mt(ntupleSvc(),"RecCgem/CgemValid");
108 StatusCode status;
109 if(cgem_ntuple){
110 status = cgem_ntuple->addItem("run", cgem_run);
111 status = cgem_ntuple->addItem("evt", cgem_evt);
112 status = cgem_ntuple->addItem("ncluster",cgem_ncluster,0,1000);
113 status = cgem_ntuple->addIndexedItem("clusterid", cgem_ncluster , cgem_clusterid);
114 status = cgem_ntuple->addIndexedItem("layerid", cgem_ncluster , cgem_layerid);
115 status = cgem_ntuple->addIndexedItem("sheetid", cgem_ncluster , cgem_sheetid);
116 status = cgem_ntuple->addIndexedItem("flag", cgem_ncluster , cgem_flag);
117 status = cgem_ntuple->addIndexedItem("energydeposit", cgem_ncluster , cgem_energydeposit);
118 status = cgem_ntuple->addIndexedItem("recphi", cgem_ncluster , cgem_recphi);
119 status = cgem_ntuple->addIndexedItem("recv", cgem_ncluster , cgem_recv);
120 status = cgem_ntuple->addIndexedItem("recZ", cgem_ncluster , cgem_recZ);
121 status = cgem_ntuple->addIndexedItem("clusterflag_first", cgem_ncluster , cgem_clusterflag_first);
122 status = cgem_ntuple->addIndexedItem("clusterflag_second", cgem_ncluster , cgem_clusterflag_second);
123 // status = cgem_ntuple->addIndexedItem("rec_ident_ID", cgem_nStrip , cgem_ident_ID);
124 status = cgem_ntuple->addItem("rec_nStrip_PerE",cgem_test_nStrip,0,3000);
125 status = cgem_ntuple->addIndexedItem("rec_ident_ID", cgem_test_nStrip , cgem_ident_ID);
126 status = cgem_ntuple->addIndexedItem("rec_strip_type", cgem_test_nStrip , cgem_strip_type);
127 status = cgem_ntuple->addIndexedItem("rec_strip_ID", cgem_test_nStrip , cgem_strip_ID);
128 status = cgem_ntuple->addIndexedItem("rec_sheet_ID", cgem_test_nStrip , cgem_sheet_ID);
129 status = cgem_ntuple->addIndexedItem("rec_layer_ID", cgem_test_nStrip , cgem_layer_ID);
130 status = cgem_ntuple->addItem("XV_ncluster",cgem_XV_ncluster,0,2000);
131 status = cgem_ntuple->addIndexedItem("rec_nStrip", cgem_XV_ncluster , cgem_nStrip);
132 status = cgem_ntuple->addIndexedItem("rec_nStrip_flagb", cgem_XV_ncluster , cgem_nStrip_flagb);
133 status = cgem_ntuple->addIndexedItem("XV_clusterid", cgem_XV_ncluster , cgem_XV_clusterid);
134 status = cgem_ntuple->addIndexedItem("XV_flag", cgem_XV_ncluster , cgem_XV_flag);
135 status = cgem_ntuple->addIndexedItem("XV_layerID", cgem_XV_ncluster , cgem_XV_layerID);
136 status = cgem_ntuple->addItem("rec_Num21",cgem_Num21,0,10000);
137 status = cgem_ntuple->addIndexedItem("rec_IsTruth21", cgem_Num21, cgem_IsTruth21);
138 status = cgem_ntuple->addIndexedItem("rec_pdg21", cgem_Num21, cgem_pdg21);
139 status = cgem_ntuple->addIndexedItem("rec_deltaZ21", cgem_Num21, cgem_deltaZ21);
140 status = cgem_ntuple->addIndexedItem("rec_deltaPhi21", cgem_Num21, cgem_deltaPhi21);
141 status = cgem_ntuple->addItem("rec_Num23",cgem_Num23,0,10000);
142 status = cgem_ntuple->addIndexedItem("rec_IsTruth23", cgem_Num23, cgem_IsTruth23);
143 status = cgem_ntuple->addIndexedItem("rec_pdg23", cgem_Num23, cgem_pdg23);
144 status = cgem_ntuple->addIndexedItem("rec_deltaZ23", cgem_Num23, cgem_deltaZ23);
145 status = cgem_ntuple->addIndexedItem("rec_deltaPhi23", cgem_Num23, cgem_deltaPhi23);
146
147 status = cgem_ntuple->addItem("MC_nhit", m_mc_nHit, 0, 1000);
148 status = cgem_ntuple->addIndexedItem("sim_nStrip",m_mc_nHit,m_mc_nStrip);
149 status = cgem_ntuple->addIndexedItem("sim_nStripB",m_mc_nHit,m_mc_nStripB);
150 status = cgem_ntuple->addIndexedItem("hittrackID",m_mc_nHit, m_mc_ID_track);
151 status = cgem_ntuple->addIndexedItem("CgemlayerID",m_mc_nHit, m_mc_ID_layer);
152 status = cgem_ntuple->addIndexedItem("pdg_code", m_mc_nHit,m_mc_pdg_code);
153 status = cgem_ntuple->addIndexedItem("ParentID", m_mc_nHit,m_mc_ID_parent);
154 status = cgem_ntuple->addIndexedItem("Deposit_Energy", m_mc_nHit,m_mc_E_deposit);
155 status = cgem_ntuple->addIndexedItem("X_pre_point", m_mc_nHit, m_mc_XYZ_pre_X );
156 status = cgem_ntuple->addIndexedItem("Y_pre_point", m_mc_nHit, m_mc_XYZ_pre_Y );
157 status = cgem_ntuple->addIndexedItem("Z_pre_point", m_mc_nHit, m_mc_XYZ_pre_Z );
158 status = cgem_ntuple->addIndexedItem("X_post_point", m_mc_nHit, m_mc_XYZ_post_X );
159 status = cgem_ntuple->addIndexedItem("Y_post_point", m_mc_nHit, m_mc_XYZ_post_Y );
160 status = cgem_ntuple->addIndexedItem("Z_post_point", m_mc_nHit, m_mc_XYZ_post_Z );
161 status = cgem_ntuple->addIndexedItem("P_X_pre_point", m_mc_nHit, m_mc_P_pre_X);
162 status = cgem_ntuple->addIndexedItem("P_Y_pre_point", m_mc_nHit, m_mc_P_pre_Y);
163 status = cgem_ntuple->addIndexedItem("P_Z_pre_point", m_mc_nHit, m_mc_P_pre_Z);
164 status = cgem_ntuple->addIndexedItem("P_X_post_point",m_mc_nHit, m_mc_P_post_X);
165 status = cgem_ntuple->addIndexedItem("P_Y_post_point",m_mc_nHit, m_mc_P_post_Y);
166 status = cgem_ntuple->addIndexedItem("P_Z_post_point",m_mc_nHit, m_mc_P_post_Z);
167 status = cgem_ntuple->addIndexedItem("rec_truth_phi", m_mc_nHit, m_mc_rec_phi);
168 status = cgem_ntuple->addIndexedItem("rec_truth_v", m_mc_nHit, m_mc_rec_v);
169 status = cgem_ntuple->addIndexedItem("rec_truth_z", m_mc_nHit, m_mc_rec_z);
170 status = cgem_ntuple->addItem("sim_nStrip_PerE",mc_test_nStrip,0,3000);
171 status = cgem_ntuple->addIndexedItem("sim_ident_ID", mc_test_nStrip , m_mc_ident_ID);
172 status = cgem_ntuple->addIndexedItem("sim_strip_type", mc_test_nStrip , m_mc_strip_type);
173 status = cgem_ntuple->addIndexedItem("sim_strip_ID" , mc_test_nStrip , m_mc_strip_ID);
174 status = cgem_ntuple->addIndexedItem("sim_sheet_ID" , mc_test_nStrip , m_mc_sheet_ID);
175 status = cgem_ntuple->addIndexedItem("sim_layer_ID" , mc_test_nStrip , m_mc_layer_ID);
176 status = cgem_ntuple->addItem("sim_Num1",m_mc_Num1,0,10000);
177 status = cgem_ntuple->addItem("sim_Num3",m_mc_Num3,0,10000);
178 status = cgem_ntuple->addItem("sim_Num21",m_mc_Num21,0,10000);
179 status = cgem_ntuple->addIndexedItem("sim_Track21", m_mc_Num21, m_mc_Track21);
180 status = cgem_ntuple->addIndexedItem("sim_PDG21", m_mc_Num21, m_mc_PDG21);
181 status = cgem_ntuple->addIndexedItem("sim_hit_PDG21", m_mc_Num21, m_mc_hit_PDG21);
182 status = cgem_ntuple->addIndexedItem("sim_deltaZ21", m_mc_Num21, m_mc_deltaZ21);
183 status = cgem_ntuple->addIndexedItem("sim_deltaPhi21", m_mc_Num21, m_mc_deltaPhi21);
184 status = cgem_ntuple->addIndexedItem("sim_hit_ID1", m_mc_Num21, m_mc_hit_ID1);
185 status = cgem_ntuple->addItem("sim_Num23",m_mc_Num23,0,10000);
186 status = cgem_ntuple->addIndexedItem("sim_Track23", m_mc_Num23, m_mc_Track23);
187 status = cgem_ntuple->addIndexedItem("sim_PDG23", m_mc_Num23, m_mc_PDG23);
188 status = cgem_ntuple->addIndexedItem("sim_hit_PDG23", m_mc_Num23, m_mc_hit_PDG23);
189 status = cgem_ntuple->addIndexedItem("sim_deltaZ23", m_mc_Num23, m_mc_deltaZ23);
190 status = cgem_ntuple->addIndexedItem("sim_deltaPhi23", m_mc_Num23, m_mc_deltaPhi23);
191 status = cgem_ntuple->addIndexedItem("sim_hit_ID2", m_mc_Num23, m_mc_hit_ID2);
192 status = cgem_ntuple->addIndexedItem("sim_hit_ID3", m_mc_Num23, m_mc_hit_ID3);
193
194 status = cgem_ntuple->addItem("MC_nTrk", m_mc_nTrk, 0, 100);
195 status = cgem_ntuple->addIndexedItem("pidTruth", m_mc_nTrk, m_mc_pidTruth);
196 status = cgem_ntuple->addIndexedItem("motheridx", m_mc_nTrk, m_mc_motheridx);
197 status = cgem_ntuple->addIndexedItem("costaTruth",m_mc_nTrk, m_mc_costaTruth);
198 status = cgem_ntuple->addIndexedItem("phiTruth", m_mc_nTrk, m_mc_phiTruth);
199 status = cgem_ntuple->addIndexedItem("Final_vrTruth", m_mc_nTrk, m_mc_Final_vrTruth);
200 status = cgem_ntuple->addIndexedItem("Final_vzTruth", m_mc_nTrk, m_mc_Final_vzTruth);
201 status = cgem_ntuple->addIndexedItem("vrTruth", m_mc_nTrk, m_mc_vrTruth);
202 status = cgem_ntuple->addIndexedItem("vzTruth", m_mc_nTrk, m_mc_vzTruth);
203 status = cgem_ntuple->addIndexedItem("pxTruth", m_mc_nTrk, m_mc_pxTruth);
204 status = cgem_ntuple->addIndexedItem("pyTruth", m_mc_nTrk, m_mc_pyTruth);
205 status = cgem_ntuple->addIndexedItem("pzTruth", m_mc_nTrk, m_mc_pzTruth);
206 status = cgem_ntuple->addIndexedItem("ptTruth", m_mc_nTrk, m_mc_ptTruth);
207 status = cgem_ntuple->addIndexedItem("pTruth", m_mc_nTrk, m_mc_pTruth);
208 status = cgem_ntuple->addIndexedItem("qTruth", m_mc_nTrk, m_mc_qTruth);
209 status = cgem_ntuple->addIndexedItem("rateTrk", m_mc_nTrk, m_mc_rateTrk);
210 status = cgem_ntuple->addIndexedItem("rateTruth", m_mc_nTrk, m_mc_rateTruth);
211 }
212
213 }
214
215 INTupleSvc* cut_ntupleSvc;
216 StatusCode src = service("NTupleSvc", cut_ntupleSvc);
217 free(foldername);
218 foldername = (char *)malloc((postfixlen + 16) * sizeof(char));
219 sprintf(foldername, "FILE145/CutResult%s", m_postfix.c_str());
220 NTuplePtr nt(ntupleSvc, foldername);
221 if ( nt ) cut_ntuple = nt;
222 else{
223 cut_ntuple = ntupleSvc->book(foldername, CLID_ColumnWiseTuple, "CgemSegment");
224 // NTuplePtr mt(ntupleSvc(),"RecCgem/CgemValid");
225 StatusCode status;
226 if(cut_ntuple){
227 status = cut_ntuple->addItem("cut_Num",cut_Num,0,10000);
228 status = cut_ntuple->addIndexedItem("cut_IsTruth21", cut_Num, cut_IsTruth21);
229 status = cut_ntuple->addIndexedItem("cut_pdg21", cut_Num, cut_pdg21);
230 status = cut_ntuple->addIndexedItem("cut_deltaZ21", cut_Num, cut_deltaZ21);
231 status = cut_ntuple->addIndexedItem("cut_deltaPhi21", cut_Num, cut_deltaPhi21);
232 status = cut_ntuple->addIndexedItem("cut_IsTruth23", cut_Num, cut_IsTruth23);
233 status = cut_ntuple->addIndexedItem("cut_pdg23", cut_Num, cut_pdg23);
234 status = cut_ntuple->addIndexedItem("cut_deltaZ23", cut_Num, cut_deltaZ23);
235 status = cut_ntuple->addIndexedItem("cut_deltaPhi23", cut_Num, cut_deltaPhi23);
236 }
237 }
238 }// if(m_checkSingleTrkMC>0)
239
240 ICgemGeomSvc* m_ICgemGeomSvc;
241 sc = service("CgemGeomSvc",m_ICgemGeomSvc);
242 if( sc != StatusCode::SUCCESS ){
243 log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
244 return StatusCode::FAILURE;
245 }
246
247
248 return StatusCode::SUCCESS;
249
250}
251
252StatusCode CgemSegmentRecAlg::execute()
253{
254 if(m_ver==1) exe_v1();
255 else if(m_ver==2) exe_v2();
256
257 return StatusCode::SUCCESS;
258}
259
260StatusCode CgemSegmentRecAlg::exe_v1()
261{
262 MsgStream log(msgSvc(), name());
263 log << MSG::INFO << "in excute()" << endreq;
264
265 //Ncut1++;
266
267 if(0 == (cgem_totevt % 10000) || m_debug>0) cout << "Event " << cgem_totevt << endl;
268 cgem_totevt++;
269
270 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
271 if (!eventHeader) {
272 log << MSG::FATAL << "Could not find Event Header" << endreq;
273 return( StatusCode::FAILURE);
274 }
275 int iEvt = eventHeader->eventNumber();
276 int iRun = eventHeader->runNumber();
277
278 bool t_ismc = false;
279 if(iRun<0){
280 t_ismc=true;
281 if(m_checkSingleTrkMC>0) getMcTruth();
282 }
283
284 double r[3] = {79.838 ,125.104 ,167.604};
285 double pi=acos(-1.0);
286 vector<double>sim_phi[3];
287 vector<double>sim_Z[3];
288 vector<int>sim_pdg[3];
289 vector<int>sim_track[3];
290 vector<int>sim_hit[3];
291
292 vector<double>phi[3];
293 vector<double>Z[3];
294 vector<int>Truth[3];
295 vector<int>rec_pdg[3];
296
297 vector<int>clusterid[3];
298 vector<int>hit_id[1000];
299 vector<int>track_id[3];
300 vector<int>trkID_S3L6;
301 // map<int,int>hit_id;
302 // map<int,int>track_id;
303 map<int,int>isATrk;
304 vector<int>pdg_track;
305 map<int,int>pdg_track_size;
306
307 m_seg_map.clear();
308 N_Seg = 0;
309
310 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
311 if (!recCgemClusterCol){
312 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
313 return StatusCode::FAILURE;
314 }
315
316 int count(0); int ss = 0;
317 int n_1DCluster=0;
318 RecCgemClusterCol::iterator itTrk=recCgemClusterCol->begin();
319 for( cgem_test_nStrip = 0, cgem_XV_ncluster = 0; itTrk != recCgemClusterCol->end(); ++itTrk)
320 {
321 if(ss>=999){N_cluster++;continue;}
322 if(cgem_test_nStrip >= 2999) {N_strip++;continue;}
323 if(cgem_XV_ncluster >= 1999) {N_XV_cluster++;continue;}
324 if((*itTrk)->getflag() == 3) N_flag3++;// combined two XV clusters
325 if(cgem_ntuple){
326
327 if((*itTrk)->getflag() == 0 || (*itTrk)->getflag() == 1) // X or V
328 {
329 int cgem_Num_Strip = (*itTrk)->getclusterflage() - (*itTrk)->getclusterflagb() + 1;
330 cgem_nStrip_flagb[cgem_XV_ncluster]=cgem_test_nStrip;
331 cgem_nStrip[cgem_XV_ncluster]=cgem_Num_Strip;
332 for(int ll = 0 ; ll < cgem_Num_Strip ; ll++){
333 cgem_ident_ID[cgem_test_nStrip] = CgemID::getIntID((*itTrk)->getlayerid(),(*itTrk)->getsheetid(),(*itTrk)->getflag(),ll + (*itTrk)->getclusterflagb());
334 cgem_layer_ID[cgem_test_nStrip] = (*itTrk)->getlayerid();
335 cgem_sheet_ID[cgem_test_nStrip] = (*itTrk)->getsheetid();
336 cgem_strip_ID[cgem_test_nStrip] = ll + (*itTrk)->getclusterflagb();
337 cgem_strip_type[cgem_test_nStrip] = (*itTrk)->getflag();
338 cgem_test_nStrip++;
339 }
340 cgem_XV_clusterid[cgem_XV_ncluster] = (*itTrk)->getclusterid();
341 cgem_XV_flag[cgem_XV_ncluster] = (*itTrk)->getflag();
342 cgem_XV_layerID[cgem_XV_ncluster] = (*itTrk)->getlayerid();
343 cgem_XV_ncluster++;
344 }
345 if((*itTrk)->getflag() == 2) // XV
346 {
347 cgem_evt = iEvt;
348 cgem_run = iRun;
349
350 cgem_clusterid[ss] = (*itTrk)->getclusterid();
351 cgem_layerid[ss] = (*itTrk)->getlayerid();
352 cgem_sheetid[ss] = (*itTrk)->getsheetid();
353 cgem_flag[ss] = (*itTrk)->getflag();
354 cgem_energydeposit[ss] = (*itTrk)->getenergydeposit();
355 cgem_recphi[ss] = (*itTrk)->getrecphi();
356 cgem_recv[ss] = (*itTrk)->getrecv();
357 cgem_recZ[ss] = (*itTrk)->getRecZ();
358 cgem_clusterflag_first[ss] = (*itTrk)->getclusterflagb();
359 cgem_clusterflag_second[ss] = (*itTrk)->getclusterflage();
360
361 }
362 }// if(cgem_ntuple)
363 if((*itTrk)->getflag() == 2)
364 {
365 phi[(*itTrk)->getlayerid()].push_back((*itTrk)->getrecphi());
366 Z[(*itTrk)->getlayerid()].push_back((*itTrk)->getRecZ());
367 ss++;
368 }
369 } // loop all clusters
370 if(cgem_ntuple) cgem_ncluster = ss;
371
372 // --- get delta_21 & delta_23
373 cgem_Num21 = 0; cgem_Num23 = 0;
374 for(int ii = 0; ii < Z[1].size(); ii++){
375 for(int jj = 0; jj < Z[0].size(); jj++){
376
377 if(cgem_Num21 == 999)continue;
378 cgem_deltaZ21[cgem_Num21] = Z[0][jj] - Z[1][ii];
379 cgem_deltaPhi21[cgem_Num21] = phi[0][jj] - phi[1][ii];
380 if(cgem_deltaPhi21[cgem_Num21]>-2*pi && cgem_deltaPhi21[cgem_Num21]<-pi)cgem_deltaPhi21[cgem_Num21] = cgem_deltaPhi21[cgem_Num21] + 2*pi;
381 if(cgem_deltaPhi21[cgem_Num21]<2*pi && cgem_deltaPhi21[cgem_Num21]>pi) cgem_deltaPhi21[cgem_Num21] = cgem_deltaPhi21[cgem_Num21] - 2*pi;
382 cgem_Num21++;
383 }
384 for(int ll = 0; ll < Z[2].size(); ll++){
385
386 if(cgem_Num23 == 999)continue;
387 cgem_deltaZ23[cgem_Num23] = Z[2][ll] - Z[1][ii];
388 cgem_deltaPhi23[cgem_Num23] = phi[2][ll] - phi[1][ii];
389 if(cgem_deltaPhi23[cgem_Num23]>-2*pi && cgem_deltaPhi23[cgem_Num23]<-pi)cgem_deltaPhi23[cgem_Num23] = cgem_deltaPhi23[cgem_Num23] + 2*pi;
390 if(cgem_deltaPhi23[cgem_Num23]<2*pi && cgem_deltaPhi23[cgem_Num23]>pi) cgem_deltaPhi23[cgem_Num23] = cgem_deltaPhi23[cgem_Num23] - 2*pi;
391 cgem_Num23++;
392 }
393 }
394
395
396 if(t_ismc){
397 MsgStream log(msgSvc(), name());
398 log << MSG::INFO << "in TruthExecute()" << endreq;
399
400 SmartDataPtr<Event::CgemMcHitCol> lv_CgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
401 if (!lv_CgemMcHitCol)
402 {
403 log << MSG::WARNING << "Could not find event" << endreq;
404 return StatusCode::FAILURE;
405 }
406
407 vector<int> hit_ID;
408 vector<int> Track_ID;
409
410 Event::CgemMcHitCol::iterator iter_truth = lv_CgemMcHitCol->begin();
411 int i=0;
412 for( mc_test_nStrip= 0; iter_truth != lv_CgemMcHitCol->end(); ++iter_truth)
413 {
414 if(i>=1000){N_hit++;continue;}
415
416 m_mc_ID_track[i]=(*iter_truth)->GetTrackID();
417 m_mc_ID_layer[i]=(*iter_truth)->GetLayerID();
418 m_mc_pdg_code[i]=(*iter_truth)->GetPDGCode();
419 m_mc_ID_parent[i]=(*iter_truth)->GetParentID();
420 m_mc_E_deposit[i]=(*iter_truth)->GetTotalEnergyDeposit();
421 m_mc_XYZ_pre_X[i] =(*iter_truth)->GetPositionXOfPrePoint() ;
422 m_mc_XYZ_pre_Y[i] =(*iter_truth)->GetPositionYOfPrePoint() ;
423 m_mc_XYZ_pre_Z[i] =(*iter_truth)->GetPositionZOfPrePoint() ;
424 m_mc_XYZ_post_X[i]=(*iter_truth)->GetPositionXOfPostPoint();
425 m_mc_XYZ_post_Y[i]=(*iter_truth)->GetPositionYOfPostPoint();
426 m_mc_XYZ_post_Z[i]=(*iter_truth)->GetPositionZOfPostPoint();
427 m_mc_P_pre_X[i]=(*iter_truth)->GetMomentumXOfPrePoint() ;
428 m_mc_P_pre_Y[i]=(*iter_truth)->GetMomentumYOfPrePoint() ;
429 m_mc_P_pre_Z[i]=(*iter_truth)->GetMomentumZOfPrePoint() ;
430 m_mc_P_post_X[i]=(*iter_truth)->GetMomentumXOfPostPoint();
431 m_mc_P_post_Y[i]=(*iter_truth)->GetMomentumYOfPostPoint();
432 m_mc_P_post_Z[i]=(*iter_truth)->GetMomentumZOfPostPoint();
433 m_mc_rec_z[i]=(m_mc_XYZ_pre_Z[i]+m_mc_XYZ_post_Z[i])/2;
434 m_mc_rec_phi[i]=(m_mc_XYZ_pre_X[i]+m_mc_XYZ_post_X[i])/2/r[m_mc_ID_layer[i]];
435
436 sim_phi[m_mc_ID_layer[i]].push_back(m_mc_rec_phi[i]);
437 sim_Z[m_mc_ID_layer[i]].push_back(m_mc_rec_z[i]);
438 sim_pdg[m_mc_ID_layer[i]].push_back(m_mc_pdg_code[i]);
439 sim_track[m_mc_ID_layer[i]].push_back(m_mc_ID_track[i]);
440 sim_hit[m_mc_ID_layer[i]].push_back(i);
441
442 ////// match to truth
443 for(int Pp =0; Pp<m_pid.size(); Pp++)
444 {
445 // if(m_mc_pdg_code[i] == m_pid[Pp] && m_mc_ID_parent[i] == 0)tmp_pdg++;
446 if(m_mc_pdg_code[i] == m_pid[Pp] ){
447 if(m_UseTesterSim){
448 if(m_mc_ID_parent[i] == 0){
449 pdg_track.push_back(m_mc_ID_track[i]);
450 pdg_track_size.insert(pair<int,int> (m_mc_ID_track[i],i));
451 }
452 }
453 else{
454 pdg_track.push_back(m_mc_ID_track[i]);
455 pdg_track_size.insert(pair<int,int> (m_mc_ID_track[i],i));
456 }
457
458 }
459 }
460 m_mc_nStripB[i] = mc_test_nStrip;
461 TArrayI ident = (*iter_truth)->GetIdentifier();
462 for(int j= 0;j<ident.GetSize();j++){
463 if(mc_test_nStrip >= 2999){N_sim_strip++; continue;}
464 m_mc_ident_ID[mc_test_nStrip] = ident.GetAt(j);
465 hit_ID.push_back(i);
466 Track_ID.push_back(m_mc_ID_track[i]);
467 mc_test_nStrip++;
468 }
469 m_mc_nStrip[i] = ident.GetSize();
470 i++;
471
472 }
473 m_mc_nHit=i;
474
475 for(m_mc_nTrk = 0; m_mc_nTrk < t_nTrkMC; m_mc_nTrk++){
476 if(m_mc_nTrk>=100){N_trk++;continue;}
477 m_mc_pidTruth[m_mc_nTrk] = t_pidTruth[m_mc_nTrk] ;
478 m_mc_motheridx[m_mc_nTrk] = t_motheridx[m_mc_nTrk] ;
479 m_mc_costaTruth[m_mc_nTrk] = t_costhetaMC[m_mc_nTrk];
480 m_mc_phiTruth[m_mc_nTrk] = t_phiMC[m_mc_nTrk] ;
481 m_mc_Final_vrTruth[m_mc_nTrk]= t_FvrMC[m_mc_nTrk] ;
482 m_mc_Final_vzTruth[m_mc_nTrk]= t_FvzMC[m_mc_nTrk] ;
483 m_mc_vrTruth[m_mc_nTrk] = t_vrMC[m_mc_nTrk] ;
484 m_mc_vzTruth[m_mc_nTrk] = t_vzMC[m_mc_nTrk] ;
485 m_mc_pxTruth[m_mc_nTrk] = t_pxMC[m_mc_nTrk] ;
486 m_mc_pyTruth[m_mc_nTrk] = t_pyMC[m_mc_nTrk] ;
487 m_mc_pzTruth[m_mc_nTrk] = t_pzMC[m_mc_nTrk] ;
488 m_mc_ptTruth[m_mc_nTrk] = t_ptMC[m_mc_nTrk] ;
489 m_mc_pTruth[m_mc_nTrk] = t_pMC[m_mc_nTrk] ;
490 m_mc_qTruth[m_mc_nTrk] = t_qMC[m_mc_nTrk] ;
491
492 for(int Pp =0; Pp<m_pid.size(); Pp++){
493 if(m_mc_pidTruth[m_mc_nTrk] == m_pid[Pp]&&fabs(m_mc_costaTruth[m_mc_nTrk])>0.93&&m_mc_ptTruth[m_mc_nTrk]<0.05)
494 OutDetecter++;
495 }
496 }
497
498 m_mc_Num21 = 0;m_mc_Num23 = 0;bool sim_IsOut21 = false;bool sim_IsOut23 = false;m_mc_Num1 = sim_Z[0].size();m_mc_Num3 = sim_Z[2].size();
499 for(int Xx = 0; Xx != sim_Z[1].size(); Xx++){
500 for(int Yy = 0; Yy != sim_Z[0].size(); Yy++){
501 if(Yy == 999){sim_IsOut21 = true;continue;}
502 m_mc_deltaZ21[m_mc_Num21] = sim_Z[0][Yy] - sim_Z[1][Xx];
503 m_mc_deltaPhi21[m_mc_Num21] = sim_phi[0][Yy] - sim_phi[1][Xx];
504 if(m_mc_deltaPhi21[m_mc_Num21]>-2*pi && m_mc_deltaPhi21[m_mc_Num21]<-pi)m_mc_deltaPhi21[m_mc_Num21] = m_mc_deltaPhi21[m_mc_Num21] + 2*pi;
505 if(m_mc_deltaPhi21[m_mc_Num21]<2*pi && m_mc_deltaPhi21[m_mc_Num21]>pi) m_mc_deltaPhi21[m_mc_Num21] = m_mc_deltaPhi21[m_mc_Num21] - 2*pi;
506
507 m_mc_hit_ID1[m_mc_Num21] = sim_hit[0][Yy];
508 m_mc_hit_ID2[m_mc_Num21] = sim_hit[1][Xx];
509 if(m_mc_pdg_code[sim_hit[1][Xx]] == m_mc_pdg_code[sim_hit[0][Yy]]) m_mc_hit_PDG21[m_mc_Num21] = m_mc_pdg_code[sim_hit[1][Xx]];
510 else m_mc_hit_PDG21[m_mc_Num21] = -999;
511 if(sim_pdg[1][Xx] == sim_pdg[0][Yy]) m_mc_PDG21[m_mc_Num21] = sim_pdg[1][Xx];
512 else m_mc_PDG21[m_mc_Num21] = -999;
513 if(sim_track[1][Xx] == sim_track[0][Yy]) m_mc_Track21[m_mc_Num21] = sim_track[1][Xx];
514 else m_mc_Track21[m_mc_Num21] = -999;
515 m_mc_Num21++;
516 }
517 for(int Zz = 0; Zz != sim_Z[2].size(); Zz++){
518 if(Zz == 999){sim_IsOut23 = true;continue;}
519 m_mc_deltaZ23[m_mc_Num23] = sim_Z[2][Zz] - sim_Z[1][Xx];
520 m_mc_deltaPhi23[m_mc_Num23] = sim_phi[2][Zz] - sim_phi[1][Xx];
521 if(m_mc_deltaPhi23[m_mc_Num23]>-2*pi && m_mc_deltaPhi23[m_mc_Num23]<-pi)m_mc_deltaPhi23[m_mc_Num23] = m_mc_deltaPhi23[m_mc_Num23] + 2*pi;
522 if(m_mc_deltaPhi23[m_mc_Num23]<2*pi && m_mc_deltaPhi23[m_mc_Num23]>pi) m_mc_deltaPhi23[m_mc_Num23] = m_mc_deltaPhi23[m_mc_Num23] - 2*pi;
523
524 m_mc_hit_ID3[m_mc_Num23] = sim_hit[2][Zz];
525 if(m_mc_pdg_code[sim_hit[1][Xx]] == m_mc_pdg_code[sim_hit[2][Zz]]) m_mc_hit_PDG23[m_mc_Num23] = m_mc_pdg_code[sim_hit[1][Xx]];
526 else m_mc_hit_PDG23[m_mc_Num23] = -999;
527 if(sim_pdg[1][Xx] == sim_pdg[2][Zz]) m_mc_PDG23[m_mc_Num23] = sim_pdg[1][Xx];
528 else m_mc_PDG23[m_mc_Num23] = -999;
529 if(sim_track[1][Xx] == sim_track[2][Zz]) m_mc_Track23[m_mc_Num23] = sim_track[1][Xx];
530 else m_mc_Track23[m_mc_Num23] = -999;
531 m_mc_Num23++;
532 }
533 }
534 if(sim_IsOut21) N_sim_N21++;
535 if(sim_IsOut23) N_sim_N23++;
536
537 for(int m = 0; m < ss; m++){
538 int IsTruth = 0;
539 map<int,int>IsEq;
540 vector<int> tmp_hit_ID[1000];
541 int tmp = 0;
542 int hit_num = 0;
543 int tmp_sim = 0;
544 // cout<<"clusterid = "<<cgem_clusterid[m]<<endl;
545
546 for(int rr=0;rr<cgem_XV_ncluster;rr++){
547 if(cgem_XV_clusterid[rr] == cgem_clusterflag_first[m] || cgem_XV_clusterid[rr] == cgem_clusterflag_second[m]){
548
549 for(int id = cgem_nStrip_flagb[rr]; id < cgem_nStrip_flagb[rr] + cgem_nStrip[rr]; id++){
550 for(int sim = 0; sim<mc_test_nStrip; sim++){
551 if(cgem_ident_ID[id]==m_mc_ident_ID[sim]){
552 tmp_hit_ID[hit_ID[sim]].push_back(Track_ID[sim]);
553 }
554 }
555 }
556 tmp += cgem_nStrip[rr];
557 }
558 }
559 for(int Mm = 0; Mm < 1000; Mm++){
560 if(tmp_hit_ID[Mm].size() == 0)continue;
561 IsEq.insert(pair<int,int> (Mm,tmp_hit_ID[Mm][0]));
562 hit_num++;
563 if(tmp_hit_ID[Mm].size() >= tmp) tmp_sim++;
564 }
565
566 bool IsPid ; int pid(0);int tmp_track_ID = 0;
567 for(map<int,int>::iterator Mm = IsEq.begin(); Mm != IsEq.end(); Mm++){
568 IsPid = false;
569 for(int Pp =0; Pp<m_pid.size(); Pp++){
570 if(m_UseTesterSim){
571 if(m_mc_pdg_code[Mm->first] == m_pid[Pp] && m_mc_ID_parent[Mm->first] == 0){
572 IsPid = true;
573 pid = m_mc_pdg_code[Mm->first];
574 }
575 }
576 else{
577 if(m_mc_pdg_code[Mm->first] == m_pid[Pp]){
578 IsPid = true;
579 pid = m_mc_pdg_code[Mm->first];
580 }
581 }
582 }
583 hit_id[m].push_back(Mm->first);
584 tmp_track_ID = Mm->second;
585 }
586 track_id[cgem_layerid[m]].push_back(tmp_track_ID);
587 clusterid[cgem_layerid[m]].push_back(m);
588 if(IsPid) rec_pdg[cgem_layerid[m]].push_back(pid);
589 else rec_pdg[cgem_layerid[m]].push_back(-999);
590 if(IsPid && tmp_sim == 1 )IsTruth =1;
591 if(IsTruth == 1)N_Match_all_hit++;
592 if(IsTruth == 1 && hit_num > 1)N_Match_hit++;
593 Truth[cgem_layerid[m]].push_back(IsTruth);
594 }
595
596 cgem_Num21 = 0;cgem_Num23 = 0;bool IsOut21 = false;bool IsOut23 = false;
597 for(int tt = 0; tt < Truth[1].size(); tt++){
598 for(int hh = 0; hh < Truth[0].size(); hh++){
599 if(cgem_Num21 == 999){IsOut21 = true;continue;}
600 if(rec_pdg[1][tt] != rec_pdg[0][hh] ) cgem_pdg21[cgem_Num21] = -9999;
601 else cgem_pdg21[cgem_Num21] = rec_pdg[1][tt];
602 if(Truth[1][tt] == 0 || Truth[0][hh] == 0)cgem_IsTruth21[cgem_Num21] = 0;
603 if(Truth[1][tt] == 1 && Truth[0][hh] == 1){
604 if(track_id[1][tt] == track_id[0][hh])cgem_IsTruth21[cgem_Num21] = 1;
605 else cgem_IsTruth21[cgem_Num21] = 0;
606 }
607 cgem_Num21++;
608 }
609
610 for(int pp = 0; pp < Truth[2].size(); pp++){
611 if(cgem_Num23 == 999){IsOut23 = true;continue;}
612 if(rec_pdg[1][tt] != rec_pdg[2][pp] ) cgem_pdg23[cgem_Num23] = -9999;
613 else cgem_pdg23[cgem_Num23] = rec_pdg[1][tt];
614 if(Truth[1][tt] == 0 || Truth[2][pp] == 0)cgem_IsTruth23[cgem_Num23] = 0;
615 if(Truth[1][tt] == 1 && Truth[2][pp] == 1){
616 if(track_id[1][tt] == track_id[2][pp])cgem_IsTruth23[cgem_Num23] = 1;
617 else cgem_IsTruth23[cgem_Num23] = 0;
618 }
619 cgem_Num23++;
620 }
621 }
622 if(IsOut21) N_N21++;
623 if(IsOut23) N_N23++;
624
625
626 int tmp_pdg[100] = {0};
627 if(pdg_track_size.size() == 0)N_PiE0++;
628 for(int Nn = 0; Nn < pdg_track.size(); Nn++){
629 int Ii = 0;
630 for(map<int,int>::iterator Mm = pdg_track_size.begin(); Mm != pdg_track_size.end(); Mm++){
631 if(pdg_track[Nn] == Mm->first) tmp_pdg[Ii]+=1;
632 Ii++;
633 }
634 }
635
636 int II = 0;
637 for( map<int,int>::iterator Mm = pdg_track_size.begin(); Mm != pdg_track_size.end(); Mm++){
638 if(tmp_pdg[II] < 3){N_PiL2++;trkID_S3L6.push_back(Mm->first);}
639 if(tmp_pdg[II] > 6){N_PiM7++;trkID_S3L6.push_back(Mm->first);}
640 II++;
641 }
642 }
643
644 if(m_UseCut){
645 N_Event++;
646
647 int cgem_Num1= phi[0].size();int cgem_Num3 = phi[2].size();
648 cut_Num = 0;
649 for(int Ii = 0; Ii < cgem_Num21; Ii++){
650 int tmp_2 = Ii/cgem_Num1;
651 int tmp_1 = Ii%cgem_Num1;
652 for(int Jj = 0; Jj < cgem_Num23; Jj++){
653 int tmp2 = Jj/cgem_Num3;
654 int tmp_3 = Jj%cgem_Num3;
655 if(tmp_2 != tmp2)continue;
656 if(!InZAr(cgem_deltaZ21[Ii],cgem_deltaZ23[Jj]) || !InPhiAr(cgem_deltaPhi21[Ii],cgem_deltaPhi23[Jj])){N_OutRange++;continue;}
657
658
659 RecCgemSegment* recsegment = new RecCgemSegment();
660 recsegment->seteventid(iEvt);
661 recsegment->setsegmentid(N_Seg);
662 recsegment->setclusterid_1(cgem_clusterid[clusterid[0][tmp_1]]);
663 recsegment->setclusterid_2(cgem_clusterid[clusterid[1][tmp_2]]);
664 recsegment->setclusterid_3(cgem_clusterid[clusterid[2][tmp_3]]);
665 int IsTruth = cgem_IsTruth21[Ii]*cgem_IsTruth23[Jj];
666 recsegment->setmatch(IsTruth);
667 m_seg_map.insert(pair<int,RecCgemSegment*>(N_Seg,recsegment));
668 N_Seg++;
669
670 cut_deltaZ21[cut_Num] = cgem_deltaZ21[Ii];
671 cut_deltaZ23[cut_Num] = cgem_deltaZ23[Jj];
672 cut_deltaPhi21[cut_Num] = cgem_deltaPhi21[Ii];
673 cut_deltaPhi23[cut_Num] = cgem_deltaPhi23[Jj];
674 cut_IsTruth21[cut_Num] = cgem_IsTruth21[Ii];
675 cut_IsTruth23[cut_Num] = cgem_IsTruth23[Jj];
676 cut_pdg21[cut_Num] = cgem_pdg21[Ii];
677 cut_pdg23[cut_Num] = cgem_pdg23[Jj];
678 cut_Num++;
679
680 if(t_ismc){
681 bool UnFit_Hits = false;
682 for(int Mm = 0; Mm <= 1000; Mm++){
683 if(Mm != clusterid[0][tmp_1] && Mm != clusterid[1][tmp_2] && Mm != clusterid[2][tmp_3]) continue;
684 for(int mm = 0; mm < hit_id[Mm].size(); mm++){
685 for(int Pp =0; Pp<trkID_S3L6.size(); Pp++){
686 if(m_mc_ID_track[hit_id[Mm][mm]] == trkID_S3L6[Pp])UnFit_Hits = true;
687 }
688 }
689 }
690 if(UnFit_Hits) continue;
691 int Un_InTrk =0;
692 if(cgem_IsTruth21[Ii]*cgem_IsTruth23[Jj] != 1){
693 int num = 0;int N_hit_size = 0;
694 for(int Mm = 0; Mm <= 1000; Mm++){
695 if(Mm != clusterid[0][tmp_1] && Mm != clusterid[1][tmp_2] && Mm != clusterid[2][tmp_3]) continue;
696 N_hit_size += hit_id[Mm].size();
697 for(int mm = 0; mm < hit_id[Mm].size(); mm++){
698 for(int Pp =0; Pp<m_pid.size(); Pp++){
699 if(m_mc_pdg_code[hit_id[Mm][mm]] == m_pid[Pp]){
700 num++;
701 if(sqrt(pow(m_mc_XYZ_pre_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_pre_Y[hit_id[Mm][mm]],2))>=sqrt(pow(m_mc_XYZ_post_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_post_Y[hit_id[Mm][mm]],2)))Un_InTrk++;
702 }
703 }
704 }
705 }
706 if(num == 0) N_nonPi++;
707 else{
708 if(num == N_hit_size){
709 if(Un_InTrk == 0 || num == Un_InTrk) N_ErrMat++;
710 else N_CrossMat++;
711 }
712 else{
713 if(num>3)N_OtherPar4++;
714 if(num==3)N_OtherPar3++;
715 if(num<3)N_OtherPar++;
716 }
717 if(num>3) N_4cluster++;
718 if(num==3) N_3cluster++;
719 }
720 continue;
721 }
722 N_Match++;
723 int InTrk = 0;
724 for(int Mm = 0; Mm <= 1000; Mm++){
725 if(Mm != clusterid[0][tmp_1] && Mm != clusterid[1][tmp_2] && Mm != clusterid[2][tmp_3]) continue;
726 for(int mm = 0; mm < hit_id[Mm].size(); mm++){
727 if(sqrt(pow(m_mc_XYZ_pre_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_pre_Y[hit_id[Mm][mm]],2))>=sqrt(pow(m_mc_XYZ_post_X[hit_id[Mm][mm]],2)+pow(m_mc_XYZ_post_Y[hit_id[Mm][mm]],2)))InTrk++;
728 }
729 }
730 if(InTrk == 3)N_inTrk3++;
731 else if(InTrk < 3&& InTrk > 0)N_InTrk++;
732 else N_outTrk++;
733 }
734 }
735 }
736
737 RecCgemSegmentCol* segmentcol = new RecCgemSegmentCol;
738
739 IDataProviderSvc* evtSvc = NULL;
740 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
741 if (evtSvc) {
742 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
743 } else {
744 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
745 return StatusCode::SUCCESS;
746 }
747
748 StatusCode segmentsc;
749 IDataManagerSvc *dataManSvc;
750 dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
751 DataObject *aSegmentCol;
752 evtSvc->findObject("/Event/Recon/RecCgemSegmentCol",aSegmentCol);
753 if(aSegmentCol != NULL) {
754 dataManSvc->clearSubTree("/Event/Recon/RecCgemSegmentCol");
755 evtSvc->unregisterObject("/Event/Recon/RecCgemSegmentCol");
756 }
757
758 segmentsc = evtSvc->registerObject("/Event/Recon/RecCgemSegmentCol", segmentcol);
759 if( segmentsc.isFailure() ) {
760 log << MSG::FATAL << "Could not register RecCgemSegment" << endreq;
761 return StatusCode::SUCCESS;
762 }
763 log << MSG::INFO << "RecCgemSegmentCol registered successfully!" <<endreq;
764
765 //push back the segment from the map//
766 for(m_seg_map_it = m_seg_map.begin();m_seg_map_it!=m_seg_map.end();++m_seg_map_it){
767 RecCgemSegment* recsegment = (*m_seg_map_it).second;
768 segmentcol->push_back(recsegment);
769 }
770
771
772 if(m_checkSingleTrkMC>0) cut_ntuple->write();
773 }
774
775 if(m_checkSingleTrkMC>0) cgem_ntuple->write();
776
777 // SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
778 // if (!recCgemSegmentCol){
779 // log << MSG::WARNING << "Could not retrieve Cgem segment list" << endreq;
780 // return StatusCode::FAILURE;
781 // }
782 //
783 // RecCgemSegmentCol::iterator itSeg = recCgemSegmentCol->begin();
784 // for(;itSeg != recCgemSegmentCol->end(); itSeg++){
785 // cout<<"event ID = "<<(*itSeg)->geteventid()<<endl;;
786 // cout<<"segment ID = "<<(*itSeg)->getsegmentid()<<endl;;
787 // cout<<"cluster ID_1 = "<<(*itSeg)->getclusterid_1()<<endl;;
788 // cout<<"cluster ID_2 = "<<(*itSeg)->getclusterid_2()<<endl;;
789 // cout<<"cluster ID_3 = "<<(*itSeg)->getclusterid_3()<<endl;;
790 // }
791
792 return StatusCode::SUCCESS;
793}
794
795StatusCode CgemSegmentRecAlg::exe_v2()
796{
797 MsgStream log(msgSvc(), name());
798
799 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
800 if (!eventHeader) {
801 log << MSG::FATAL << "Could not find Event Header" << endreq;
802 return( StatusCode::FAILURE);
803 }
804 m_run = eventHeader->runNumber();
805 m_evt = eventHeader->eventNumber();
806 cgem_totevt++;
807
808 resetData();
809
810 /* --- get cluster information */
811 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
812 if (!recCgemClusterCol){
813 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
814 return StatusCode::FAILURE;
815 }
816
817
818 if(m_debug) {
819 cout<<"get "<<recCgemClusterCol->size()<<" Cgem clusters:"<<endl;
820 cout<<setw(15)<<"idClu"<<setw(20)<<"phi"<<setw(20)<<"Z"<<endl;
821 }
822 RecCgemClusterCol::iterator itCluster=recCgemClusterCol->begin();
823 for( ; itCluster != recCgemClusterCol->end(); ++itCluster)
824 {
825 int CluFlag = (*itCluster)->getflag();
826 if(CluFlag==2)
827 {
828 int idClu = (*itCluster)->getclusterid();
829 int iLayer = (*itCluster)->getlayerid();
830 int iSheet = (*itCluster)->getsheetid();
831 int id_XClu = (*itCluster)->getclusterflagb();
832 int id_VClu = (*itCluster)->getclusterflage();
833 double recPhi = (*itCluster)->getrecphi();
834 double recV = (*itCluster)->getrecv();
835 double recZ = (*itCluster)->getRecZ();
836 mVec_idXVClu[iLayer].push_back(idClu);
837 mVec_phiClu[iLayer].push_back(recPhi);
838 mVec_idXClu[iLayer].push_back(id_XClu);
839 mVec_ZClu[iLayer].push_back(recZ);
840 mVec_idVClu[iLayer].push_back(id_VClu);
841 if(m_debug) cout<<setw(15)<<idClu<<setw(20)<<recPhi<<setw(20)<<recZ<<endl;
842 }
843 }
844
845 for(int i=0; i<3; i++) n_XVCluster[i] = mVec_phiClu[i].size();
846 /* get cluster information --- */
847
848
849 // --- register RecCgemSegmentCol
850 RecCgemSegmentCol* segmentcol = new RecCgemSegmentCol;
851 IDataProviderSvc* evtSvc = NULL;
852 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
853 if (evtSvc) {
854 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
855 } else {
856 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
857 return StatusCode::SUCCESS;
858 }
859
860 StatusCode segmentsc;
861 IDataManagerSvc *dataManSvc;
862 dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
863 DataObject *aSegmentCol;
864 evtSvc->findObject("/Event/Recon/RecCgemSegmentCol",aSegmentCol);
865 if(aSegmentCol != NULL) {
866 dataManSvc->clearSubTree("/Event/Recon/RecCgemSegmentCol");
867 evtSvc->unregisterObject("/Event/Recon/RecCgemSegmentCol");
868 }
869
870 segmentsc = evtSvc->registerObject("/Event/Recon/RecCgemSegmentCol", segmentcol);
871 if( segmentsc.isFailure() ) {
872 log << MSG::FATAL << "Could not register RecCgemSegment" << endreq;
873 return StatusCode::SUCCESS;
874 }
875 log << MSG::INFO << "RecCgemSegmentCol registered successfully!" <<endreq;
876
877
878 /* --- search good segment candidates */
879 int i_GoodSeg = 0;
880 for(int i=0; i<n_XVCluster[1]; i++) // layer 2
881 {
882 int nGoodSegCandi=0;
883 double phi2 = mVec_phiClu[1][i];
884 double Z2 = mVec_ZClu[1][i];
885 int id_lay2 = mVec_idXVClu[1][i];
886 vector<double> vec_dPhi21;
887 vector<double> vec_dPhi23;
888 vector<double> vec_dZ21;
889 vector<double> vec_dZ23;
890 for(int j=0; j<n_XVCluster[0]; j++) // layer 1
891 {
892 double phi1 = mVec_phiClu[0][j];
893 double Z1 = mVec_ZClu[0][j];
894 double dPhi21 = phi1-phi2;
895 while(dPhi21<-(CLHEP::pi)) dPhi21+=CLHEP::twopi;
896 while(dPhi21> (CLHEP::pi)) dPhi21-=CLHEP::twopi;
897 double dZ21 = Z1-Z2;
898 vec_dPhi21.push_back(dPhi21);
899 vec_dZ21.push_back(dZ21);
900 }
901 for(int k=0; k<n_XVCluster[2]; k++) // layer 3
902 {
903 int id_lay3 = mVec_idXVClu[2][k];
904 double phi3 = mVec_phiClu[2][k];
905 double Z3 = mVec_ZClu[2][k];
906 double dPhi23 = phi3-phi2;
907 while(dPhi23<-(CLHEP::pi)) dPhi23+=CLHEP::twopi;
908 while(dPhi23> (CLHEP::pi)) dPhi23-=CLHEP::twopi;
909 double dZ23 = Z3-Z2;
910 vec_dPhi23.push_back(dPhi23);
911 vec_dZ23.push_back(dZ23);
912 for(int j=0; j<n_XVCluster[0]; j++)
913 {
914 int id_lay1 = mVec_idXVClu[0][j];
915 double dPhi21 = vec_dPhi21[j];
916 double dZ21 = vec_dZ21[j];
917 bool isGoodDPhi = InPhiAr(dPhi21,dPhi23);
918 bool isGoodDZ = InZAr(dZ21,dZ23);
919 if(isGoodDPhi && isGoodDZ) {
920 nGoodSegCandi++;
921 RecCgemSegment* recsegment = new RecCgemSegment();
922 recsegment->seteventid(m_evt);
923 recsegment->setsegmentid(i_GoodSeg);
924 recsegment->setclusterid_1(id_lay1);
925 recsegment->setclusterid_2(id_lay2);
926 recsegment->setclusterid_3(id_lay3);
927 //int IsTruth = cgem_IsTruth21[Ii]*cgem_IsTruth23[Jj];
928 //recsegment->setmatch(IsTruth);
929 segmentcol->push_back(recsegment);
930 if(m_debug) {
931 cout<<"id_lay1, id_lay2, id_lay3 = "<<id_lay1<<", "<<id_lay2<<", "<<id_lay3<<endl;
932 cout<<"getclusterid_3: "<<recsegment->getclusterid_3()<<endl;
933 }
934 i_GoodSeg++;
935 }
936 }
937 }
938
939 if(m_debug) cout<<"nGoodSegCandi="<<nGoodSegCandi<<" for XV cluster "<<i<<" in layer 2"<<endl;
940 }// loop XV-clusters on layer 2
941 /* --- search good segment candidates */
942
943 if(m_debug) printInfo();
944
945 return StatusCode::SUCCESS;
946}
947
948void CgemSegmentRecAlg::resetData()
949{
950 for(int i=0; i<3; i++)
951 {
952 mVec_idXVClu[i].clear();
953 mVec_phiClu[i].clear();
954 mVec_idXClu[i].clear();
955 mVec_ZClu[i].clear();
956 mVec_idVClu[i].clear();
957 //mVec_dPhi21.clear();
958 //mVec_dPhi23.clear();
959 //mVec_dZ21.clear();
960 //mVec_dZ23.clear();
961 }
962}
963
964void CgemSegmentRecAlg::printInfo()
965{
966 cout<<"CgemSegmentRecAlg printInfo(): "<<endl;
967 cout.precision(4);
968 int txtWidth=10;
969 for(int i=0; i<3; i++)
970 {
971 int nXV = mVec_phiClu[i].size();
972 cout<<"Layer "<<i<<endl;
973 cout<<setw(txtWidth)<<"phi"<<setw(txtWidth)<<"Z"<<setw(txtWidth)<<"id_X"<<setw(txtWidth)<<"id_V"<<endl;
974 for(int j=0; j<nXV; j++)
975 {
976 cout<<setw(txtWidth)<<mVec_phiClu[i][j]<<setw(txtWidth)<<mVec_ZClu[i][j]<<setw(txtWidth)<<mVec_idXClu[i][j]<<setw(txtWidth)<<mVec_idVClu[i][j]<<endl;
977 }
978 }
979
980 SmartDataPtr<RecCgemSegmentCol> recCgemSegmentCol(eventSvc(), "/Event/Recon/RecCgemSegmentCol");
981 if(recCgemSegmentCol)
982 {
983 int n_seg = recCgemSegmentCol->size();
984 RecCgemSegmentCol::iterator it_seg = recCgemSegmentCol->begin();
985 cout<<"RecCgemSegmentCol contains "<<n_seg<<" CgemSegments: "<<endl;
986 cout<<setw(txtWidth)<<"seg_id"<<setw(txtWidth)<<"clu_id1"<<setw(txtWidth)<<"clu_id2"<<setw(txtWidth)<<"clu_id3"<<endl;
987 for(; it_seg!=recCgemSegmentCol->end(); it_seg++)
988 {
989 //RecCgemSegment* recsegment = (RecCgemSegment*) (*it_seg);
990 //DstCgemSegment* recsegment = (RecCgemSegment*) (*it_seg);
991 cout <<setw(txtWidth)<<(*it_seg)->getsegmentid()
992 <<setw(txtWidth)<<(*it_seg)->getclusterid_1()
993 <<setw(txtWidth)<<(*it_seg)->getclusterid_2()
994 <<setw(txtWidth)<<(*it_seg)->getclusterid_3()
995 //<<setw(txtWidth)<<(*it_seg)->geteventid()
996 <<endl;
997 //cout<< recsegment->getcluserid_3()<<endl;
998 }
999 }
1000
1001 cout.precision(6);
1002}
1003
1004void CgemSegmentRecAlg::getMcTruth(){
1005 MsgStream log(msgSvc(), name());
1006 StatusCode sc;
1007 t_t0Truth=-1;
1008 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
1009 if (!mcParticleCol) {
1010 log << MSG::ERROR<< "Could not find McParticle" << endreq;
1011 } else{
1012 int rootIndex(-1);
1013 int itk= 0;
1014 t_nTrkMC=0;
1015 bool findJpsi=false;
1016 if(m_debug>0) cout<< "size_mcParticleCol "<<mcParticleCol->size()<< endl;
1017 McParticleCol::iterator iter_mc = mcParticleCol->begin();
1018 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
1019 //if(!(*iter_mc)->primaryParticle()) continue;
1020 int pdg=(*iter_mc)->particleProperty();
1021 if(pdg==443){
1022 findJpsi=true;
1023 rootIndex = (*iter_mc)->trackIndex();
1024 }
1025 // if(!findJpsi) continue;
1026
1027 t_pidTruth[itk] = (*iter_mc)->particleProperty();
1028 if(m_debug>0) cout<< "tk "<<itk<< " pid="<< t_pidTruth[itk]<< endl;
1029
1030 int mcidx = ((*iter_mc)->mother()).trackIndex() - rootIndex;
1031 if (mcidx > 0) t_motheridx[itk] = mcidx-1;
1032 else t_motheridx[itk] = mcidx;
1033
1034
1035 //if((m_pid!=0) && (abs(t_pidTruth[itk]) != m_pid)){ continue; }
1036
1037 t_t0Truth=(*iter_mc)->initialPosition().t();
1038
1039 if(itk>=100) break;
1040
1041 int pid = t_pidTruth[itk];
1042 if( pid == 0 ) {
1043 log << MSG::WARNING << "Wrong particle id" <<endreq;
1044 t_qMC[itk]=0;
1045 }else{
1046 IPartPropSvc* p_PartPropSvc;
1047 static const bool CREATEIFNOTTHERE(true);
1048 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
1049 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
1050 cout<< "CgemSegmentRecAlg::getMcTruth() Could not initialize Particle Properties Service" << endl;
1051 }
1052 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
1053 std::string name;
1054 if( p_particleTable->particle(pid) ){
1055 name = p_particleTable->particle(pid)->name();
1056 t_qMC[itk] = p_particleTable->particle(pid)->charge();
1057 }else if( p_particleTable->particle(-pid) ){
1058 name = "anti " + p_particleTable->particle(-pid)->name();
1059 t_qMC[itk] = (-1.)*p_particleTable->particle(-pid)->charge();
1060 }
1061 }
1062
1063 t_vrMC[itk] = sqrt((*iter_mc)->initialPosition().x()*(*iter_mc)->initialPosition().x() + (*iter_mc)->initialPosition().y()*(*iter_mc)->initialPosition().y()) ;
1064 t_vzMC[itk] = (*iter_mc)->initialPosition().z();
1065 t_pxMC[itk] = (*iter_mc)->initialFourMomentum().px();
1066 t_pyMC[itk] = (*iter_mc)->initialFourMomentum().py();
1067 t_pzMC[itk] = (*iter_mc)->initialFourMomentum().pz();
1068 t_ptMC[itk] = sqrt(t_pxMC[itk]*t_pxMC[itk] + t_pyMC[itk]*t_pyMC[itk]);
1069 t_pMC[itk] = sqrt(t_ptMC[itk]*t_ptMC[itk] + t_pzMC[itk]*t_pzMC[itk]);
1070 t_FvrMC[itk]= sqrt((*iter_mc)->finalPosition().x()*(*iter_mc)->finalPosition().x() + (*iter_mc)->finalPosition().y()*(*iter_mc)->finalPosition().y()) ;
1071 t_FvzMC[itk] = (*iter_mc)->initialPosition().z();
1072
1073
1074 Hep3Vector p3(t_pxMC[itk],t_pyMC[itk],t_pzMC[itk]);
1075
1076 t_costhetaMC[itk] = t_pzMC[itk]/t_pMC[itk];
1077
1078 t_phiMC[itk] = p3.phi();
1079 /*
1080 m_mc_p = t_pMC[itk];
1081 m_mc_pt = t_ptMC[itk];
1082 m_mc_vr = t_vrMC[itk];
1083 m_mc_vz = t_vzMC[itk];
1084 m_mc_phi = t_phiMC[itk];
1085 m_mc_ct = t_costhetaMC[itk];
1086 m_ttuple->write();
1087
1088 if(m_debug>1){
1089 cout<<"Truth tk "<<itk <<" pid "<<pid
1090 <<" pxyz("<< t_pxMC[itk]<<"," <<t_pyMC[itk]<<","<< t_pzMC[itk]
1091 <<") pt " <<t_ptMC[itk]<<" p "<< t_pMC[itk]
1092 <<" cosTa "<<t_costhetaMC[itk]
1093 <<" phi "<<t_phiMC[itk]<<" charge "<<t_qMC[itk] <<endl;
1094 }
1095 */
1096 itk++;
1097 t_nTrkMC++;
1098 }
1099 }
1100
1101 // //mc hit
1102 // SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
1103 // if (!mdcMcHitCol) {
1104 // log << MSG::ERROR<< "Could not find MdcMcHitCol" << endreq;
1105 // } else{
1106 // int ihit= 0;
1107 // t_nMdcMcHit = mdcMcHitCol->size();
1108 // if(m_debug>1) cout<< "nMcHit="<<t_nMdcMcHit << endl;
1109 // Event::MdcMcHitCol::iterator iter_mchit = mdcMcHitCol->begin();
1110 // for (;iter_mchit != mdcMcHitCol->end(); iter_mchit++,ihit++ ) {
1111 // const Identifier id= (*iter_mchit)->identify();
1112 // //int layer = MdcID::layer(id);
1113 // //int wire = MdcID::wire(id);
1114 // int iMcTk = (*iter_mchit)->getTrackIndex();
1115 // double mcX = (*iter_mchit)->getPositionX()/10.;
1116 // double mcY = 3*iter_mchit)->getPositionY()/10.;
1117 // double mcZ = (*iter_mchit)->getPositionZ()/10.;
1118 // double mcLR = (*iter_mchit)->getPositionFlag();
1119 // double mcDrift = (*iter_mchit)->getDriftDistance()/10.;
1120 // double mcEnergy = (*iter_mchit)->getDepositEnergy()/10.;
1121 // if (mcLR == 0) mcLR = -1;
1122 // if(m_debug){
1123 // cout<<"Truth hit "<<ihit <<" tk "<< iMcTk
1124 // <<" xyz("<< mcX<<"," <<mcY<<","<< mcZ
1125 // <<") lr " <<mcLR
1126 // <<" drift " <<mcDrift
1127 // <<" energy " <<mcEnergy
1128 // <<endl;
1129 // }
1130 // }
1131 // }
1132
1133}
1134
1135bool CgemSegmentRecAlg::InZAr(double x, double y) {
1136 if(fabs(x)>m_z_bound||fabs(y)>m_z_bound)return false;
1137
1138 double distance,ZCut,deltaL;
1139 bool IsIn = false;
1140 //if(m_debug>0) cout<<"m_z_slope.size() = "<<m_z_slope.size()<<endl;
1141 for(int z = 0; z < m_z_slope.size(); z++){
1142 distance = (y - m_z_slope[z]*x)/sqrt(1 + pow(m_z_slope[z],2));
1143 deltaL = (x + m_z_slope[z]*y)/sqrt(1 + pow(m_z_slope[z],2));
1144 ZCut = m_z_nSigma*(m_z_Par1[z]*pow(deltaL,2) + m_z_Par2[z]*deltaL + m_z_Par3[z]);
1145 if(fabs(distance) < ZCut) IsIn = true;
1146 }
1147 return IsIn;
1148}
1149
1150bool CgemSegmentRecAlg::InPhiAr(double x, double y) {
1151 if(fabs(x)>m_phi_bound||fabs(y)>m_phi_bound)return false;
1152
1153 bool IsIn = false;
1154 double distance = (y - m_phi_slope*x)/sqrt(1 + pow(m_phi_slope,2));
1155 double deltaL = (x + m_phi_slope*y)/sqrt(1 + pow(m_phi_slope,2));
1156 double mean = m_phi_3Par[0]*pow(deltaL,3) + m_phi_3Par[1]*pow(deltaL,2) + m_phi_3Par[2]*deltaL + m_phi_3Par[3];
1157 double PhiCut = m_phi_nSigma*(m_phi_4Par[0]*pow(deltaL,4) + m_phi_4Par[1]*pow(deltaL,3) + m_phi_4Par[2]*pow(deltaL,2) + m_phi_4Par[3]*deltaL + m_phi_4Par[4]);
1158 if(fabs(distance - mean) < PhiCut) IsIn = true;
1159 return IsIn;
1160}
1161
1162
1163// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1164StatusCode CgemSegmentRecAlg::finalize() {
1165
1166 MsgStream log(msgSvc(), name());
1167 log << MSG::INFO << "in finalize()" << endreq;
1168
1169 //TFile fout(m_fname.c_str(), m_fmode.c_str());
1170 //fout.Close();
1171
1172 if(m_debug>0)
1173 {
1174 cout << "Cgem total event: " << cgem_totevt << endl;
1175
1176 cout << "Out of Detector track: " << OutDetecter << endl;
1177 cout << "Number of hits which IsTruth equare to 1: " << N_Match_all_hit << endl;
1178 cout << "Number of hits which IsTruth equare to 1 and matched hits more than 2: " << N_Match_hit << endl;
1179
1180 cout << "Cgem events of ncluster more than 1000 : "<<N_cluster<<endl;
1181 cout << "Cgem events of XV_ncluster more than 2000 : "<<N_XV_cluster<<endl;
1182 cout << "Cgem events of reconstruction strips more than 3000 : "<<N_strip<<endl;
1183 cout << "Cgem events of cluster with falg = 3 : "<<N_flag3<<endl;
1184 cout << "Cgem events of hits more than 1000 : "<<N_hit<<endl;
1185 cout << "Cgem events of simulation strips more than 3000 : "<<N_sim_strip<<endl;
1186 cout << "Cgem events of tracks more than 100 : "<<N_trk<<endl;
1187 cout << "Cgem events of rec Num_21 more than 1000 : "<<N_N21<<endl;
1188 cout << "Cgem events of rec Num_23 more than 1000 : "<<N_N23<<endl;
1189 cout << "Cgem events of sim Num_21 more than 1000 : "<<N_sim_N21<<endl;
1190 cout << "Cgem events of sim Num_23 more than 1000 : "<<N_sim_N23<<endl;
1191 cout << "Cgem segments stay out of range : "<<N_OutRange<<endl;
1192 cout << "Cgem events of particles = 0 : "<<N_PiE0<<endl;
1193 cout << "Cgem tracks of particles between 0 to 3 : "<<N_PiL2<<endl;
1194 cout << "Cgem tracks of particles more than 6 : "<<N_PiM7<<endl;
1195 cout << "Cgem events' number : "<<N_Event<<endl;
1196 cout << "Cgem segments' number : "<<N_Seg<<endl;
1197 cout << "Matching Cgem segments' number : "<<N_Match<<endl;
1198 cout << "Matching Cgem segments' number which towards out of detector : "<<N_outTrk<<endl;
1199 cout << "Matching Cgem segments' number which towards in detector : "<<N_inTrk3<<endl;
1200 cout << "Matching Cgem segments' number which contains other particles towards in detector : "<<N_InTrk<<endl;
1201 cout << "Unmatched Cgem segments' number with equare to 3 pid clusters : "<<N_3cluster<<endl;
1202 cout << "Unmatched Cgem segments' number with more than 3 pid clusters : "<<N_4cluster<<endl;
1203 cout << "Unmatched Cgem segments' number without pid clusters : "<<N_nonPi<<endl;
1204 cout << "Unmatched Cgem segments' number which reconstructed error : "<<N_ErrMat<<endl;
1205 cout << "Unmatched Cgem segments' number which both contains toward in and out track : "<<N_CrossMat<<endl;
1206 cout << "Unmatched Cgem segments' number which contains more than 3 pid clusters effects by other particles : "<<N_OtherPar4<<endl;
1207 cout << "Unmatched Cgem segments' number which contains equare to 3 pid clusters effects by other particles : "<<N_OtherPar3<<endl;
1208 cout << "Unmatched Cgem segments' number which contains 1 or 2 pid clusters and effects by other particles : "<<N_OtherPar<<endl;
1209
1210
1211 cout<<endl<<"******* CgemSegmentRecAlg *******"<<endl;
1212 cout<<"* m_dst ="<<m_dst<<endl;
1213 cout<<"* nEvent ="<<Ncut0<<endl;
1214 cout<<"***************************"<<endl;
1215 cout<<"******** continueHit ******"<<endl;
1216
1217 // for(int i=0;i<43;i++){
1218 // cout<<"* lay"<<i<<" size"<<m_continuousHit[i].size()<<" element:";
1219 // for(int j=0;j<continuousHit[i].size();j++)
1220 // cout<<" "<<continuousHit[i][j];
1221 // cout<<endl;
1222 // }
1223
1224 }
1225
1226 return StatusCode::SUCCESS;
1227}
phi2
Double_t phi2
Definition: DataBase/tau_mode.c:7
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
phi1
Double_t phi1
Definition: DataBase/tau_mode.c:7
RecCgemSegmentCol
ObjectVector< RecCgemSegment > RecCgemSegmentCol
Definition: RecCgemSegment.h:44
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
CgemID::getIntID
static value_type getIntID(unsigned int f_layer, unsigned int f_sheet, unsigned int f_strip_type, unsigned int f_strip)
Definition: CgemID.cxx:104
CgemSegmentRecAlg::CgemSegmentRecAlg
CgemSegmentRecAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition: CgemSegmentRecAlg.cxx:55
CgemSegmentRecAlg::InZAr
bool InZAr(double x, double y)
Definition: CgemSegmentRecAlg.cxx:1135
CgemSegmentRecAlg::InPhiAr
bool InPhiAr(double x, double y)
Definition: CgemSegmentRecAlg.cxx:1150
DstCgemSegment::getclusterid_3
const int getclusterid_3() const
Definition: DstCgemSegment.h:38
DstCgemSegment::setclusterid_2
void setclusterid_2(const int clusterId_2)
Definition: DstCgemSegment.h:48
DstCgemSegment::setmatch
void setmatch(const int match)
Definition: DstCgemSegment.h:52
DstCgemSegment::setsegmentid
void setsegmentid(const int segmentId)
Definition: DstCgemSegment.h:50
DstCgemSegment::setclusterid_1
void setclusterid_1(const int clusterId_1)
Definition: DstCgemSegment.h:47
DstCgemSegment::seteventid
void seteventid(const int eventId)
Definition: DstCgemSegment.h:51
DstCgemSegment::setclusterid_3
void setclusterid_3(const int clusterId_3)
Definition: DstCgemSegment.h:49
Event
Definition: Event.h:21
num
int num[96]
Definition: ranlxd.c:373
pi
const float pi
Definition: vector3.h:133