12#include "GaudiKernel/AlgFactory.h"
13#include "GaudiKernel/DataObject.h"
14#include "GaudiKernel/IEventProcessor.h"
16#include "GaudiKernel/Incident.h"
17#include "GaudiKernel/IIncidentSvc.h"
18#include "GaudiKernel/Memory.h"
23#include "GaudiKernel/ISvcLocator.h"
24#include "GaudiKernel/IDataProviderSvc.h"
25#include "GaudiKernel/Bootstrap.h"
26#include "GaudiKernel/RegistryEntry.h"
27#include "GaudiKernel/MsgStream.h"
38#include "GaudiKernel/SmartDataPtr.h"
47 Algorithm(name,pSvcLocator){
49 declareProperty(
"selGoodDigi",m_selGoodDigi=1);
50 declareProperty(
"minDigiTime",m_minDigiTime=-8875);
51 declareProperty(
"maxDigiTime",m_maxDigiTime=-8562);
52 declareProperty(
"minQDigi",myQMin=0.0);
53 declareProperty(
"LUTfile", lutfile =
"/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_from_17_to_17.root");
54 declareProperty(
"align_flag", align_flag =
false);
69void TestTrack::DefineHitTree() {
71 tree->Branch(
"event", &event,
"event/I");
72 tree->Branch(
"nhit", &nhit,
"nhit/I");
74 tree->Branch(
"hit_strip", &hit_strip,
"hit_strip[nhit]/I");
75 tree->Branch(
"hit_view", &hit_view,
"hit_view[nhit]/I");
76 tree->Branch(
"hit_layer", &hit_layer,
"hit_layer[nhit]/I");
77 tree->Branch(
"hit_sheet", &hit_sheet,
"hit_sheet[nhit]/I");
78 tree->Branch(
"hit_length", &hit_length,
"hit_length[nhit]/D");
80 tree->Branch(
"hit_channel", &hit_channel,
"hit_channel[nhit]/I");
81 tree->Branch(
"hit_roc", &hit_roc,
"hit_roc[nhit]/I");
82 tree->Branch(
"hit_feb", &hit_feb,
"hit_feb[nhit]/I");
83 tree->Branch(
"hit_tiger", &hit_tiger,
"hit_tiger[nhit]/I");
84 tree->Branch(
"hit_chip", &hit_chip,
"hit_chip[nhit]/I");
86 tree->Branch(
"hit_t", &hit_t,
"hit_t[nhit]/D");
87 tree->Branch(
"hit_q", &hit_q,
"hit_q[nhit]/D");
88 tree->Branch(
"hit_quality", &hit_quality,
"hit_quality[nhit]/I");
89 tree->Branch(
"hit_selgooddigi", &hit_selgooddigi,
"hit_selgooddigi[nhit]/I");
94void TestTrack::DefineClusterTree() {
96 tree->Branch(
"ncluster", &ncluster,
"ncluster/I");
98 tree->Branch(
"ncluster_1d", &ncluster_1d,
"ncluster_1d/I");
99 tree->Branch(
"ncluster_1d_L1_S1_x", &ncluster_1d_L1_S1_x,
"ncluster_1d_L1_S1_x/I");
100 tree->Branch(
"ncluster_1d_L2_S1_x", &ncluster_1d_L2_S1_x,
"ncluster_1d_L2_S1_x/I");
101 tree->Branch(
"ncluster_1d_L2_S2_x", &ncluster_1d_L2_S2_x,
"ncluster_1d_L2_S2_x/I");
102 tree->Branch(
"ncluster_1d_L1_S1_v", &ncluster_1d_L1_S1_v,
"ncluster_1d_L1_S1_v/I");
103 tree->Branch(
"ncluster_1d_L2_S1_v", &ncluster_1d_L2_S1_v,
"ncluster_1d_L2_S1_v/I");
104 tree->Branch(
"ncluster_1d_L2_S2_v", &ncluster_1d_L2_S2_v,
"ncluster_1d_L2_S2_v/I");
105 tree->Branch(
"cluster_1d_ID", &cluster_1d_ID,
"cluster_1d_ID[ncluster]/D");
106 tree->Branch(
"cluster_1d_t", &cluster_1d_t,
"cluster_1d_t[ncluster]/D");
107 tree->Branch(
"cluster_1d_q", &cluster_1d_q,
"cluster_1d_q[ncluster]/D");
108 tree->Branch(
"cluster_1d_r", &cluster_1d_r,
"cluster_1d_r[ncluster]/D");
109 tree->Branch(
"cluster_1d_v", &cluster_1d_v,
"cluster_1d_v[ncluster]/D");
110 tree->Branch(
"cluster_1d_v_cc", &cluster_1d_v_cc,
"cluster_1d_v_cc[ncluster]/D");
111 tree->Branch(
"cluster_1d_v_tpc", &cluster_1d_v_tpc,
"cluster_1d_v_tpc[ncluster]/D");
112 tree->Branch(
"cluster_1d_phi", &cluster_1d_phi,
"cluster_1d_phi[ncluster]/D");
113 tree->Branch(
"cluster_1d_phi_cc", &cluster_1d_phi_cc,
"cluster_1d_phi_cc[ncluster]/D");
114 tree->Branch(
"cluster_1d_phi_tpc", &cluster_1d_phi_tpc,
"cluster_1d_phi_tpc[ncluster]/D");
115 tree->Branch(
"cluster_1d_layerid", &cluster_1d_layerid,
"cluster_1d_layerid[ncluster]/I");
116 tree->Branch(
"cluster_1d_sheetid", &cluster_1d_sheetid,
"cluster_1d_sheetid[ncluster]/I");
117 tree->Branch(
"cluster_1d_view", &cluster_1d_view,
"cluster_1d_view[ncluster]/I");
118 tree->Branch(
"cluster_1d_strip1", &cluster_1d_strip1,
"cluster_1d_strip1[ncluster]/I");
119 tree->Branch(
"cluster_1d_strip2", &cluster_1d_strip2,
"cluster_1d_strip2[ncluster]/I");
120 tree->Branch(
"cluster_1d_size", &cluster_1d_size,
"cluster_1d_size[ncluster]/I");
121 tree->Branch(
"cluster_1d_a_tpc", &cluster_1d_a_tpc,
"cluster_1d_a_tpc[ncluster]/D");
122 tree->Branch(
"cluster_1d_b_tpc", &cluster_1d_b_tpc,
"cluster_1d_b_tpc[ncluster]/D");
123 tree->Branch(
"cluster_1d_hitindex", cluster_1d_hitindex,
"cluster_1d_hitindex[ncluster][50]/I");
126 tree->Branch(
"ncluster_2d", &ncluster_2d,
"ncluster_2d/I");
127 tree->Branch(
"cluster_2d_t", &cluster_2d_t,
"cluster_2d_t[ncluster]/D");
128 tree->Branch(
"cluster_2d_q", &cluster_2d_q,
"cluster_2d_q[ncluster]/D");
129 tree->Branch(
"cluster_2d_r", &cluster_2d_r,
"cluster_2d_r[ncluster]/D");
130 tree->Branch(
"cluster_2d_z", &cluster_2d_z,
"cluster_2d_z[ncluster]/D");
131 tree->Branch(
"cluster_2d_z_cc", &cluster_2d_z_cc,
"cluster_2d_z_cc[ncluster]/D");
132 tree->Branch(
"cluster_2d_z_tpc", &cluster_2d_z_tpc,
"cluster_2d_z_tpc[ncluster]/D");
133 tree->Branch(
"cluster_2d_phi", &cluster_2d_phi,
"cluster_2d_phi[ncluster]/D");
134 tree->Branch(
"cluster_2d_phi_cc", &cluster_2d_phi_cc,
"cluster_2d_phi_cc[ncluster]/D");
135 tree->Branch(
"cluster_2d_phi_tpc", &cluster_2d_phi_tpc,
"cluster_2d_phi_tpc[ncluster]/D");
136 tree->Branch(
"cluster_2d_layerid", &cluster_2d_layerid,
"cluster_2d_layerid[ncluster]/I");
137 tree->Branch(
"cluster_2d_sheetid", &cluster_2d_sheetid,
"cluster_2d_sheetid[ncluster]/I");
138 tree->Branch(
"cluster_2d_view", &cluster_2d_view,
"cluster_2d_view[ncluster]/I");
139 tree->Branch(
"cluster_2d_idx", &cluster_2d_idx,
"cluster_2d_idx[ncluster]/I");
140 tree->Branch(
"cluster_2d_idv", &cluster_2d_idv,
"cluster_2d_idv[ncluster]/I");
141 tree->Branch(
"cluster_2d_highest", &cluster_2d_highest,
"cluster_2d_highest[ncluster]/I");
146void TestTrack::DefineTrackTree() {
149 tree->Branch(
"track_dr", &track_dr,
"track_dr/D");
150 tree->Branch(
"track_phi0", &track_phi0,
"track_phi0/D");
151 tree->Branch(
"track_dz", &track_dz,
"track_dz/D");
152 tree->Branch(
"track_tanL", &track_tanL,
"track_tanL/D");
153 tree->Branch(
"track_chi2", &track_chi2,
"track_chi2/D");
156 tree->Branch(
"track_nfitpoint", &track_nfitpoint,
"track_nfitpoint/I");
157 tree->Branch(
"track_clusterid", &track_clusterid,
"track_clusterid[track_nfitpoint]/I");
158 tree->Branch(
"track_layerid", &track_layerid,
"track_layerid[track_nfitpoint]/I");
159 tree->Branch(
"track_sheetid", &track_sheetid,
"track_sheetid[track_nfitpoint]/I");
162 tree->Branch(
"track_test_layerid", &track_test_layerid,
"track_test_layerid/I");
163 tree->Branch(
"track_test_sheetid", &track_test_sheetid,
"track_test_sheetid/I");
166 tree->Branch(
"track_xpoca_glo", &track_xpoca_glo,
"track_xpoca_glo/D");
167 tree->Branch(
"track_ypoca_glo", &track_ypoca_glo,
"track_ypoca_glo/D");
168 tree->Branch(
"track_zpoca_glo", &track_zpoca_glo,
"track_zpoca_glo/D");
172 tree->Branch(
"track_x1top_glo", &track_x1top_glo,
"track_x1top_glo/D");
173 tree->Branch(
"track_y1top_glo", &track_y1top_glo,
"track_y1top_glo/D");
174 tree->Branch(
"track_z1top_glo", &track_z1top_glo,
"track_z1top_glo/D");
175 tree->Branch(
"track_phi1top_loc", &track_phi1top_loc,
"track_phi1top_loc/D");
176 tree->Branch(
"track_v1top_loc", &track_v1top_loc,
"track_v1top_loc/D");
177 tree->Branch(
"track_x1bot_glo", &track_x1bot_glo,
"track_x1bot_glo/D");
178 tree->Branch(
"track_y1bot_glo", &track_y1bot_glo,
"track_y1bot_glo/D");
179 tree->Branch(
"track_z1bot_glo", &track_z1bot_glo,
"track_z1bot_glo/D");
180 tree->Branch(
"track_phi1bot_loc", &track_phi1bot_loc,
"track_phi1bot_loc/D");
181 tree->Branch(
"track_v1bot_loc", &track_v1bot_loc,
"track_v1bot_loc/D");
183 tree->Branch(
"track_x2top_glo", &track_x2top_glo,
"track_x2top_glo/D");
184 tree->Branch(
"track_y2top_glo", &track_y2top_glo,
"track_y2top_glo/D");
185 tree->Branch(
"track_z2top_glo", &track_z2top_glo,
"track_z2top_glo/D");
186 tree->Branch(
"track_phi2top_loc", &track_phi2top_loc,
"track_phi2top_loc/D");
187 tree->Branch(
"track_v2top_loc", &track_v2top_loc,
"track_v2top_loc/D");
188 tree->Branch(
"track_x2bot_glo", &track_x2bot_glo,
"track_x2bot_glo/D");
189 tree->Branch(
"track_y2bot_glo", &track_y2bot_glo,
"track_y2bot_glo/D");
190 tree->Branch(
"track_z2bot_glo", &track_z2bot_glo,
"track_z2bot_glo/D");
191 tree->Branch(
"track_phi2bot_loc", &track_phi2bot_loc,
"track_phi2bot_loc/D");
192 tree->Branch(
"track_v2bot_loc", &track_v2bot_loc,
"track_v2bot_loc/D");
195 tree->Branch(
"ang_xy_L1", &ang_xy_L1,
"ang_xy_L1/D");
196 tree->Branch(
"ang_yz_L1", &ang_yz_L1,
"ang_yz_L1/D");
198 tree->Branch(
"ang_xy_L2", &ang_xy_L2,
"ang_xy_L2/D");
199 tree->Branch(
"ang_yz_L2", &ang_yz_L2,
"ang_yz_L2/D");
218 for(
int ihit = 0; ihit <
MAXNOFHITS; ihit++) {
219 hit_strip[ihit] = -1;
221 hit_layer[ihit] = -1;
222 hit_sheet[ihit] = -1;
223 hit_length[ihit] = -1;
224 hit_channel[ihit] = -1;
227 hit_tiger[ihit] = -1;
231 hit_quality[ihit] = 0;
232 hit_selgooddigi[ihit] = 0;
239 ncluster_1d_L1_S1_x = 0;
240 ncluster_1d_L2_S1_x = 0;
241 ncluster_1d_L2_S2_x = 0;
242 ncluster_1d_L1_S1_v = 0;
243 ncluster_1d_L2_S1_v = 0;
244 ncluster_1d_L2_S2_v = 0;
248 cluster_1d_t[iclu] = 0.;
249 cluster_1d_q[iclu] = 0.;
250 cluster_1d_r[iclu] = 0.;
251 cluster_1d_v[iclu] = 0.;
252 cluster_1d_phi[iclu] = 0.;
253 cluster_1d_v_cc[iclu] = 0.;
254 cluster_1d_phi_cc[iclu] = 0.;
255 cluster_1d_v_tpc[iclu] = 0.;
256 cluster_1d_phi_tpc[iclu] = 0.;
257 cluster_1d_a_tpc[iclu] = 0.;
258 cluster_1d_b_tpc[iclu] = 0.;
259 cluster_1d_layerid[iclu] = -1;
260 cluster_1d_sheetid[iclu] = -1;
261 cluster_1d_view[iclu] = -1;
262 cluster_1d_size[iclu] = -1;
263 cluster_1d_strip1[iclu] = -1;
264 cluster_1d_strip2[iclu] = -1;
265 for(
int ihit = 0; ihit <
MAXNOFHITS; ihit++) cluster_1d_hitindex[iclu][ihit] = -1;
275 ncluster_2d_L1_S1 = 0;
276 ncluster_2d_L2_S1 = 0;
277 ncluster_2d_L2_S2 = 0;
281 cluster_2d_t[iclu] = 0.;
282 cluster_2d_q[iclu] = 0.;
283 cluster_2d_r[iclu] = 0.;
284 cluster_2d_z[iclu] = 0.;
285 cluster_2d_phi[iclu] = 0.;
286 cluster_2d_z_cc[iclu] = 0.;
287 cluster_2d_phi_cc[iclu] = 0.;
288 cluster_2d_z_tpc[iclu] = 0.;
289 cluster_2d_phi_tpc[iclu] = 0.;
290 cluster_2d_layerid[iclu] = -1;
291 cluster_2d_sheetid[iclu] = -1;
292 cluster_2d_view[iclu] = -1;
293 cluster_2d_idx[iclu] = -1;
294 cluster_2d_idv[iclu] = -1;
295 cluster_2d_highest[iclu] = -1;
300 map_L1_S1_stripx_to_hit.clear();
301 map_L2_S1_stripx_to_hit.clear();
302 map_L2_S2_stripx_to_hit.clear();
303 map_L1_S1_stripv_to_hit.clear();
304 map_L2_S1_stripv_to_hit.clear();
305 map_L2_S2_stripv_to_hit.clear();
320 track_clusterid[iclus] = -1;
321 track_layerid[iclus] = -1;
322 track_sheetid[iclus] = -1;
326 track_test_layerid = -1;
327 track_test_sheetid = -1;
330 track_xpoca_glo = -9999;
331 track_ypoca_glo = -9999;
332 track_zpoca_glo = -9999;
336 track_x1top_glo = -9999;
337 track_y1top_glo = -9999;
338 track_z1top_glo = -9999;
339 track_x1bot_glo = -9999;
340 track_y1bot_glo = -9999;
341 track_z1bot_glo = -9999;
343 track_x2top_glo = -9999;
344 track_y2top_glo = -9999;
345 track_z2top_glo = -9999;
346 track_x2bot_glo = -9999;
347 track_y2bot_glo = -9999;
348 track_z2bot_glo = -9999;
353bool TestTrack::selDigi(CgemDigiCol::iterator
iter,
int sel)
355 if(sel==0)
return true;
357 double time = (*iter)->getTime_ns();
358 bool timeIsGood=
true;
359 if(time<m_minDigiTime||time>m_maxDigiTime) timeIsGood=
false;
361 double Q = (*iter)->getCharge_fc();
362 bool chargeIsGood=
true;
363 if(Q<myQMin) chargeIsGood=
false;
371 if(is_xstrip ==
true) view = 0;
372 int quality = lutreader->
GetQuality(layer, sheet, view, strip);
373 bool qualityIsGood=
true;
374 if(quality!=0) qualityIsGood=
false;
376 if(sel==1)
return timeIsGood&&chargeIsGood;
377 else if(sel==-1)
return !timeIsGood&&chargeIsGood;
378 else if(sel==2)
return timeIsGood&&chargeIsGood&&qualityIsGood;
385std::map<int, std::vector < int > > * TestTrack::GetStripToHitMap(
int ilayer,
int isheet,
int iview) {
388 if(ilayer == 0 && isheet == 0)
return &map_L1_S1_stripx_to_hit;
389 else if(ilayer == 1 && isheet == 0)
return &map_L2_S2_stripx_to_hit;
390 else if(ilayer == 1 && isheet == 1)
return &map_L2_S1_stripx_to_hit;
393 if(ilayer == 0 && isheet == 0)
return &map_L1_S1_stripv_to_hit;
394 else if(ilayer == 1 && isheet == 0)
return &map_L2_S2_stripv_to_hit;
395 else if(ilayer == 1 && isheet == 1)
return &map_L2_S1_stripv_to_hit;
403 MsgStream log(
msgSvc(), name());
404 log << MSG::INFO <<
"TestTrack initialize()" << endreq;
409 sc = service(
"CgemGeomSvc", m_SvcCgem);
410 if(sc != StatusCode::SUCCESS) {
411 log << MSG::ERROR <<
"can not use CgemGeomSvc" << endreq;
412 return StatusCode::FAILURE;
428 output =
new TFile(
"track.root",
"RECREATE");
429 tree =
new TTree(
"tree",
"hit, cluster and track info tree");
436 return StatusCode::SUCCESS;
442 MsgStream log(
msgSvc(), name());
448 ISvcLocator* svcLocator = Gaudi::svcLocator();
449 StatusCode sc=svcLocator->service(
"EventDataSvc", m_evtSvc);
451 cout<<
"Could not accesss EventDataSvc!"<<endl;
454 SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,
"/Event/Digi/CgemDigiCol");
456 cout<<
"Could not retrieve CGEM digi collection"<<endl;
457 nhit = aDigiCol->size();
463 CgemDigiCol::iterator iDigiCol;
464 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
466 const Identifier ident = (*iDigiCol)->identify();
472 if(is_xstrip ==
true) view = 0;
473 double charge = (*iDigiCol)->getCharge_fc();
474 double time = (*iDigiCol)->getTime_ns();
476 hit_strip[ihit] = strip;
477 hit_layer[ihit] = layer;
478 hit_sheet[ihit] = sheet;
479 hit_view[ihit] = view;
481 hit_q[ihit] = charge;
485 else hit_length[ihit] = anode->
getLength();
487 int channel = lutreader->
GetChannel(layer, sheet, view, strip);
488 int roc = lutreader->
GetROC(layer, sheet, view, strip);
489 int feb = lutreader->
GetFEB(layer, sheet, view, strip);
490 int tiger = lutreader->
GetTIGER(layer, sheet, view, strip);
491 int chip = lutreader->
GetChip(layer, sheet, view, strip);
492 int quality = lutreader->
GetQuality(layer, sheet, view, strip);
494 hit_channel[ihit] = channel;
497 hit_tiger[ihit] = tiger;
498 hit_chip[ihit] = chip;
499 hit_quality[ihit] = quality;
501 if(selDigi(iDigiCol, m_selGoodDigi)==
true) {
502 hit_selgooddigi[ihit] = 1;
506 std::map<int, std::vector < int > > *map_strip_to_hit = GetStripToHitMap(layer, sheet, view);
507 std::map<int, std::vector < int > >::iterator it1 = map_strip_to_hit->find(strip);
508 if (it1 != map_strip_to_hit->end()) {
509 it1->second.push_back(ihit);
512 std::vector< int > hitlist;
513 hitlist.push_back(ihit);
514 std::pair<int, std::vector < int> > pair(strip, hitlist);
515 map_strip_to_hit->insert(pair);
526 SmartDataPtr<RecCgemClusterCol> aClusterCol(m_evtSvc,
"/Event/Recon/RecCgemClusterCol");
528 cout<<
"Could not retrieve CGEM cluster collection"<<endl;
529 int nclu = aClusterCol->size();
532 std::cout <<
"nclu " << nclu << std::endl;
540 int tmp_cluster_L1_S1 = -1;
541 int tmp_cluster_L2_S1 = -1;
542 int tmp_cluster_L2_S2 = -1;
543 double tmp_charge_L1_S1 = 0.;
544 double tmp_charge_L2_S1 = 0.;
545 double tmp_charge_L2_S2 = 0.;
549 RecCgemClusterCol::iterator iClusterCol;
550 for(iClusterCol=aClusterCol->begin(); iClusterCol!=aClusterCol->end(); iClusterCol++)
555 int flag = (*iClusterCol)->getflag();
556 int clusterid = (*iClusterCol)->getclusterid();
557 int trkid = (*iClusterCol)->getTrkId();
558 int layerid = (*iClusterCol)->getlayerid();
559 int sheetid = (*iClusterCol)->getsheetid();
560 double edep = (*iClusterCol)->getenergydeposit();
561 double phi = (*iClusterCol)->getrecphi();
562 double v = (*iClusterCol)->getrecv();
563 double z = (*iClusterCol)->getRecZ();
564 double phi_cc = (*iClusterCol)->getrecphi_CC();
565 double v_cc = (*iClusterCol)->getrecv_CC();
566 double z_cc = (*iClusterCol)->getRecZ_CC();
567 double phi_tpc = (*iClusterCol)->getrecphi_mTPC();
568 double v_tpc = (*iClusterCol)->getrecv_mTPC();
569 double z_tpc = (*iClusterCol)->getRecZ_mTPC();
585 int flagB = (*iClusterCol)->getclusterflagb();
586 int flagE = (*iClusterCol)->getclusterflage();
589 if(flag==0 || flag == 1) {
590 cluster_1d_ID[ncluster] = clusterid;
592 cluster_1d_t[ncluster] = 0;
593 cluster_1d_q[ncluster] = edep;
595 cluster_1d_layerid[ncluster] = layerid;
596 cluster_1d_sheetid[ncluster] = sheetid;
597 cluster_1d_view[ncluster] = flag;
598 cluster_1d_r[ncluster] = anode_mid_gap;
599 cluster_1d_phi[ncluster] = phi;
600 cluster_1d_phi_cc[ncluster] = phi_cc;
601 cluster_1d_phi_tpc[ncluster] = phi_tpc;
602 cluster_1d_v[ncluster] =
v;
603 cluster_1d_v_cc[ncluster] = v_cc;
604 cluster_1d_v_tpc[ncluster] = v_tpc;
605 cluster_1d_a_tpc[ncluster] = a_tpc;
606 cluster_1d_b_tpc[ncluster] = b_tpc;
608 cluster_1d_strip1[ncluster] = flagB;
609 cluster_1d_strip2[ncluster] = flagE;
610 int size = flagE - flagB + 1;
611 cluster_1d_size[ncluster] = size;
614 std::map<int, std::vector < int > > *read_map_strip_to_hit = GetStripToHitMap(layerid, sheetid, flag);
615 std::map<int, std::vector < int > >::iterator it2;
624 for(
int istrip = 0; istrip < size; istrip++) {
625 int stripid = flagB + istrip;
626 it2 = read_map_strip_to_hit->find(stripid);
627 std::vector< int > hitlist = it2->second;
628 int size_hitlist = hitlist.size();
629 cluster_1d_hitindex[ncluster][istrip] = hitlist[size_hitlist-1];
635 else if(flag==2 || flag==3) {
636 cluster_2d_ID[ncluster] = clusterid;
638 cluster_2d_t[ncluster] = 0;
639 cluster_2d_q[ncluster] = edep;
641 cluster_2d_layerid[ncluster] = layerid;
642 cluster_2d_sheetid[ncluster] = sheetid;
643 if(layerid==0 && phi>0) cluster_2d_sheetid[ncluster] = 1;
644 cluster_2d_view[ncluster] = flag;
645 cluster_2d_r[ncluster] = anode_mid_gap;
646 cluster_2d_z[ncluster] = z;
647 cluster_2d_z_cc[ncluster] = z_cc;
648 cluster_2d_z_tpc[ncluster] = z_tpc;
649 cluster_2d_phi[ncluster] = phi;
650 cluster_2d_phi_cc[ncluster] = phi_cc;
651 cluster_2d_phi_tpc[ncluster] = phi_tpc;
653 cluster_2d_idx[ncluster] = flagB;
654 cluster_2d_idv[ncluster] = flagE;
665 if(edep > tmp_charge_L1_S1) {
666 tmp_charge_L1_S1 = edep; tmp_cluster_L1_S1 = ncluster;
669 else if(layerid==1 && sheetid==0) {
670 if(edep > tmp_charge_L2_S1){
671 tmp_charge_L2_S1 = edep; tmp_cluster_L2_S1 = ncluster;
674 else if(layerid==1 && sheetid==1) {
675 if(edep > tmp_charge_L2_S2){
676 tmp_charge_L2_S2 = edep; tmp_cluster_L2_S2 = ncluster;
687 if(layerid == 0) ncluster_1d_L1_S1_x++;
688 else if(layerid == 1) {
689 if(sheetid == 0) ncluster_1d_L2_S1_x++;
690 else ncluster_1d_L2_S2_x++;
695 if(layerid == 0) ncluster_1d_L1_S1_v++;
696 else if(layerid == 1) {
697 if(sheetid == 0) ncluster_1d_L2_S1_v++;
698 else ncluster_1d_L2_S2_v++;
703 if(layerid == 0) ncluster_2d_L1_S1++;
704 else if(layerid == 1) {
705 if(sheetid == 0) ncluster_2d_L2_S1++;
706 else ncluster_2d_L2_S2++;
716 if(tmp_cluster_L1_S1!=-1) cluster_2d_highest[tmp_cluster_L1_S1]=1;
717 if(tmp_cluster_L2_S1!=-1) cluster_2d_highest[tmp_cluster_L2_S1]=1;
718 if(tmp_cluster_L2_S2!=-1) cluster_2d_highest[tmp_cluster_L2_S2]=1;
722 SmartDataPtr<RecMdcTrackCol> aTrackCol(m_evtSvc,
"/Event/Recon/RecMdcTrackCol");
729 RecMdcTrackCol::iterator iTrackCol;
730 for(iTrackCol=aTrackCol->begin(); iTrackCol!=aTrackCol->end(); iTrackCol++)
742 if(track->
trackId() != 0)
continue;
744 track_chi2 = track->
chi2();
746 track_dr = track->
helix(0);
747 track_phi0 = track->
helix(1);
748 track_dz = track->
helix(3);
749 track_tanL = track->
helix(4);
752 track_nfitpoint = vecclusters.size();
755 bool istestplane[2][2];
756 for(
int ilay=0; ilay<
MAXNOFLAYER; ilay++)
for(
int ishe=0; ishe<
MAXNOFSHEET; ishe++) istestplane[ilay][ishe] =
true;
758 for(
int iclus=0; iclus < track_nfitpoint; iclus++) {
761 double z = tcluster->
getRecZ();
764 track_layerid[iclus] = tcluster->
getlayerid();
765 track_sheetid[iclus] = tcluster->
getsheetid();
766 if(phi > 0) track_sheetid[iclus] = 1;
771 istestplane[track_layerid[iclus]][track_sheetid[iclus]] =
false;
779 if(istestplane[ilay][ishe] ==
true) {
780 track_test_layerid = ilay;
781 track_test_sheetid = ishe;
789 double xP;
double yP;
double zP;
790 ComputePOCA(xP, yP, zP);
793 double xP1;
double yP1;
double zP1;
double phiP1;
double vP1;
794 double xP2;
double yP2;
double zP2;
double phiP2;
double vP2;
796 double angCR_xy_L1, angCR_yz_L1;
798 bool gotit1 = ComputeIntersection(0, xP1, yP1, zP1, phiP1, vP1, xP2, yP2, zP2, phiP2, vP2, angCR_xy_L1, angCR_yz_L1);
800 ang_xy_L1 = angCR_xy_L1;
801 ang_yz_L1 = angCR_yz_L1;
804 double xP3;
double yP3;
double zP3;
double phiP3;
double vP3;
805 double xP4;
double yP4;
double zP4;
double phiP4;
double vP4;
807 double angCR_xy_L2, angCR_yz_L2;
809 bool gotit2 = ComputeIntersection(1, xP3, yP3, zP3, phiP3, vP3, xP4, yP4, zP4, phiP4, vP4, angCR_xy_L2, angCR_yz_L2);
811 ang_xy_L2 = angCR_xy_L2;
812 ang_yz_L2 = angCR_yz_L2;
814 if(gotit1==
false || gotit2==
false) cout <<
"TestTrack: error in intersection calculation. Intersection on L1 " << gotit1 <<
", on L2 " << gotit2 << endl;
818 track_xpoca_glo = xP; track_ypoca_glo = yP; track_zpoca_glo = zP;
821 track_x1top_glo = xP1; track_y1top_glo = yP1; track_z1top_glo = zP1; track_phi1top_loc = phiP1; track_v1top_loc = vP1;
822 track_x1bot_glo = xP2; track_y1bot_glo = yP2; track_z1bot_glo = zP2; track_phi1bot_loc = phiP2; track_v1bot_loc = vP2;
824 track_x2top_glo = xP3; track_y2top_glo = yP3; track_z2top_glo = zP3; track_phi2top_loc = phiP3; track_v2top_loc = vP3;
825 track_x2bot_glo = xP4; track_y2bot_glo = yP4; track_z2bot_glo = zP4; track_phi2bot_loc = phiP4; track_v2bot_loc = vP4;
834 return StatusCode::SUCCESS;
838 MsgStream log(
msgSvc(),name());
839 log << MSG::INFO <<
"TestTrack finalize()" << endreq;
842 return StatusCode::SUCCESS;
848void TestTrack::ComputePOCA(
double &xP,
double &yP,
double &zP)
851 xP = track_dr *
cos(track_phi0);
852 yP = track_dr *
sin(track_phi0);
861bool TestTrack::ComputeIntersection(
int layerid,
862 double &x1,
double &y1,
double &z1,
double&
phi1,
double &v1,
863 double &x2,
double &y2,
double &z2,
double&
phi2,
double &v2,
864 double &angCR_xy,
double &angCR_yz)
866 StraightLine linefit(track_dr, track_phi0, track_dz, track_tanL);
867 HepPoint3D posup, posdown, posup_temp, posdown_temp;
868 double phivup[2], phivdown[2], phivup_temp[2], phivdown_temp[2];
871 if(align_flag==
true) {
873 gotit = midplane->
getPointAligned(layerid, 0, linefit, posup_temp, posdown, phivup_temp, phivdown);
875 gotit = midplane->
getPointAligned(layerid, 1, linefit, posup, posdown_temp, phivup, phivdown_temp);
877 else gotit = midplane->
getPointIdealGeom(layerid, linefit, posup, posdown, phivup, phivdown);
891 Hep3Vector CosmicRay_xy(x1-x2, y1-y2, 0);
892 Hep3Vector CosmicRay_yz( 0, y1-y2, z1-z2);
894 Hep3Vector NormalToCRayInPoint(x1, y1, 0);
895 Hep3Vector zAxis(0, 0, 1);
897 angCR_xy = CosmicRay_xy.angle(NormalToCRayInPoint);
898 angCR_yz = CosmicRay_yz.angle(zAxis);
double sin(const BesAngle a)
double cos(const BesAngle a)
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
SmartRefVector< RecCgemCluster > ClusterRefVec
double getVStripLength(int V_ID) const
static int strip(const Identifier &id)
static int sheet(const Identifier &id)
static int layer(const Identifier &id)
static bool is_xstrip(const Identifier &id)
int GetChip(int ilayer, int isheet, int iview, int istrip)
int GetQuality(int ilayer, int isheet, int iview, int istrip)
int GetROC(int ilayer, int isheet, int iview, int istrip)
int GetChannel(int ilayer, int isheet, int iview, int istrip)
int GetTIGER(int ilayer, int isheet, int iview, int istrip)
int GetFEB(int ilayer, int isheet, int iview, int istrip)
bool getPointIdealGeom(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
bool getPointAligned(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
const double chi2() const
const int trackId() const
const HepVector helix() const
......
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
virtual CgemGeoAlign * getAlignPtr() const =0
virtual CgemMidDriftPlane * getMidDriftPtr() const =0
double getRecZ(void) const
int getclusterid(void) const
int getlayerid(void) const
double getrecphi(void) const
int getsheetid(void) const
const ClusterRefVec getVecClusters() const
TestTrack(const std::string &name, ISvcLocator *pSvcLocator)