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"
37#include "GaudiKernel/SmartDataPtr.h"
45 Algorithm(name,pSvcLocator){
63 MsgStream log(
msgSvc(), name());
64 log << MSG::INFO <<
"TestCluster initialize()" << endreq;
69 sc = service(
"CgemGeomSvc", m_SvcCgem);
70 if(sc != StatusCode::SUCCESS) {
71 log << MSG::ERROR <<
"can not use CgemGeomSvc" << endreq;
72 return StatusCode::FAILURE;
75 output =
new TFile(
"cluster.root",
"RECREATE");
76 tree =
new TTree(
"tree",
"cluster info tree");
78 tree->Branch(
"ncluster", &ncluster,
"ncluster/I");
79 tree->Branch(
"ncluster_L1_S1_x", &ncluster_L1_S1_x,
"ncluster_L1_S1_x/I");
80 tree->Branch(
"ncluster_L2_S1_x", &ncluster_L2_S1_x,
"ncluster_L2_S1_x/I");
81 tree->Branch(
"ncluster_L2_S2_x", &ncluster_L2_S2_x,
"ncluster_L2_S2_x/I");
82 tree->Branch(
"ncluster_L1_S1_v", &ncluster_L1_S1_v,
"ncluster_L1_S1_v/I");
83 tree->Branch(
"ncluster_L2_S1_v", &ncluster_L2_S1_v,
"ncluster_L2_S1_v/I");
84 tree->Branch(
"ncluster_L2_S2_v", &ncluster_L2_S2_v,
"ncluster_L2_S2_v/I");
85 tree->Branch(
"anode_radius_L1_x", &anode_radius_L1_x,
"anode_radius_L1_x/D");
86 tree->Branch(
"anode_radius_L1_v", &anode_radius_L1_v,
"anode_radius_L1_v/D");
87 tree->Branch(
"anode_mid_gap_L1", &anode_mid_gap_L1,
"anode_mid_gap_L1/D");
88 tree->Branch(
"anode_radius_L2_x", &anode_radius_L2_x,
"anode_radius_L2_x/D");
89 tree->Branch(
"anode_radius_L2_v", &anode_radius_L2_v,
"anode_radius_L2_v/D");
90 tree->Branch(
"anode_mid_gap_L2", &anode_mid_gap_L2,
"anode_mid_gap_L2/D");
91 tree->Branch(
"cluster_t", &cluster_t,
"cluster_t[ncluster]/D");
92 tree->Branch(
"cluster_q", &cluster_q,
"cluster_q[ncluster]/D");
93 tree->Branch(
"cluster_qx", &cluster_qx,
"cluster_qx[ncluster]/D");
94 tree->Branch(
"cluster_qv", &cluster_qv,
"cluster_qv[ncluster]/D");
95 tree->Branch(
"cluster_r", &cluster_r,
"cluster_r[ncluster]/D");
96 tree->Branch(
"cluster_z", &cluster_z,
"cluster_z[ncluster]/D");
97 tree->Branch(
"cluster_z_cc", &cluster_z_cc,
"cluster_z_cc[ncluster]/D");
98 tree->Branch(
"cluster_z_tpc", &cluster_z_tpc,
"cluster_z_tpc[ncluster]/D");
99 tree->Branch(
"cluster_v", &cluster_v,
"cluster_v[ncluster]/D");
100 tree->Branch(
"cluster_v_cc", &cluster_v_cc,
"cluster_v_cc[ncluster]/D");
101 tree->Branch(
"cluster_v_tpc", &cluster_v_tpc,
"cluster_v_tpc[ncluster]/D");
102 tree->Branch(
"cluster_phi", &cluster_phi,
"cluster_phi[ncluster]/D");
103 tree->Branch(
"cluster_phi_cc", &cluster_phi_cc,
"cluster_phi_cc[ncluster]/D");
104 tree->Branch(
"cluster_phi_tpc", &cluster_phi_tpc,
"cluster_phi_tpc[ncluster]/D");
105 tree->Branch(
"cluster_layerid", &cluster_layerid,
"cluster_layerid[ncluster]/I");
106 tree->Branch(
"cluster_sheetid", &cluster_sheetid,
"cluster_sheetid[ncluster]/I");
107 tree->Branch(
"cluster_view", &cluster_view,
"cluster_view[ncluster]/I");
108 tree->Branch(
"cluster_idx", &cluster_idx,
"cluster_idx[ncluster]/I");
109 tree->Branch(
"cluster_idv", &cluster_idv,
"cluster_idv[ncluster]/I");
110 tree->Branch(
"cluster_stripx1", &cluster_stripx1,
"cluster_stripx1[ncluster]/I");
111 tree->Branch(
"cluster_stripx2", &cluster_stripx2,
"cluster_stripx2[ncluster]/I");
112 tree->Branch(
"cluster_stripv1", &cluster_stripv1,
"cluster_stripv1[ncluster]/I");
113 tree->Branch(
"cluster_stripv2", &cluster_stripv2,
"cluster_stripv2[ncluster]/I");
114 tree->Branch(
"cluster_sizex", &cluster_sizex,
"cluster_sizex[ncluster]/I");
115 tree->Branch(
"cluster_sizev", &cluster_sizev,
"cluster_sizev[ncluster]/I");
116 tree->Branch(
"cluster_highest", &cluster_highest,
"cluster_highest[ncluster]/I");
117 tree->Branch(
"cluster_ax_tpc", &cluster_ax_tpc,
"cluster_ax_tpc[ncluster]/D");
118 tree->Branch(
"cluster_bx_tpc", &cluster_bx_tpc,
"cluster_bx_tpc[ncluster]/D");
119 tree->Branch(
"cluster_av_tpc", &cluster_av_tpc,
"cluster_av_tpc[ncluster]/D");
120 tree->Branch(
"cluster_bv_tpc", &cluster_bv_tpc,
"cluster_bv_tpc[ncluster]/D");
122 return StatusCode::SUCCESS;
127 ncluster_L1_S1_x = 0;
128 ncluster_L2_S1_x = 0;
129 ncluster_L2_S2_x = 0;
130 ncluster_L1_S1_v = 0;
131 ncluster_L2_S1_v = 0;
132 ncluster_L2_S2_v = 0;
136 cluster_t[iclu] = 0.;
137 cluster_q[iclu] = 0.;
138 cluster_qx[iclu] = 0.;
139 cluster_qv[iclu] = 0.;
140 cluster_r[iclu] = 0.;
141 cluster_z[iclu] = 0.;
142 cluster_v[iclu] = 0.;
143 cluster_phi[iclu] = 0.;
144 cluster_z_cc[iclu] = 0.;
145 cluster_v_cc[iclu] = 0.;
146 cluster_phi_cc[iclu] = 0.;
147 cluster_z_tpc[iclu] = 0.;
148 cluster_v_tpc[iclu] = 0.;
149 cluster_phi_tpc[iclu] = 0.;
150 cluster_ax_tpc[iclu] = 0.;
151 cluster_bx_tpc[iclu] = 0.;
152 cluster_av_tpc[iclu] = 0.;
153 cluster_bv_tpc[iclu] = 0.;
154 cluster_layerid[iclu] = 0;
155 cluster_sheetid[iclu] = 0;
156 cluster_view[iclu] = 0;
157 cluster_sizex[iclu] = 0;
158 cluster_sizev[iclu] = 0;
159 cluster_idx[iclu] = 0;
160 cluster_idv[iclu] = 0;
161 cluster_stripx1[iclu] = 0;
162 cluster_stripx2[iclu] = 0;
163 cluster_stripv1[iclu] = 0;
164 cluster_stripv2[iclu] = 0;
165 cluster_highest[iclu] = 0;
171 MsgStream log(
msgSvc(), name());
175 ISvcLocator* svcLocator = Gaudi::svcLocator();
176 StatusCode sc=svcLocator->service(
"EventDataSvc", m_evtSvc);
178 cout<<
"Could not accesss EventDataSvc!"<<endl;
181 SmartDataPtr<RecCgemClusterCol> aClusterCol(m_evtSvc,
"/Event/Recon/RecCgemClusterCol");
183 cout<<
"Could not retrieve CGEM cluster collection"<<endl;
185 int nclu = aClusterCol->size();
188 std::cout <<
"nclu " << nclu << std::endl;
199 anode_radius_L1_x = anode->
getRX();
200 anode_radius_L1_v = anode->
getRV();
203 anode_radius_L2_x = anode->
getRX();
204 anode_radius_L2_v = anode->
getRV();
208 int tmp_cluster_L1_S1 = -1;
209 int tmp_cluster_L2_S1 = -1;
210 int tmp_cluster_L2_S2 = -1;
211 double tmp_charge_L1_S1 = 0.;
212 double tmp_charge_L2_S1 = 0.;
213 double tmp_charge_L2_S2 = 0.;
215 RecCgemClusterCol::iterator iClusterCol;
216 for(iClusterCol=aClusterCol->begin(); iClusterCol!=aClusterCol->end(); iClusterCol++)
221 int flag = (*iClusterCol)->getflag();
224 if(flag==2 || flag==3) nclusterxv++;
226 int clusterid = (*iClusterCol)->getclusterid();
227 int trkid = (*iClusterCol)->getTrkId();
228 int layerid = (*iClusterCol)->getlayerid();
229 int sheetid = (*iClusterCol)->getsheetid();
231 double edep = (*iClusterCol)->getenergydeposit();
232 double phi = (*iClusterCol)->getrecphi();
233 double v = (*iClusterCol)->getrecv();
234 double z = (*iClusterCol)->getRecZ();
235 double phi_cc = (*iClusterCol)->getrecphi_CC();
236 double v_cc = (*iClusterCol)->getrecv_CC();
237 double z_cc = (*iClusterCol)->getRecZ_CC();
238 double phi_tpc = (*iClusterCol)->getrecphi_mTPC();
239 double v_tpc = (*iClusterCol)->getrecv_mTPC();
240 double z_tpc = (*iClusterCol)->getRecZ_mTPC();
251 cluster_t[ncluster] = 0;
253 cluster_q[ncluster] = edep;
257 double anode_radius_x = anode->
getRX();
258 double anode_radius_v = anode->
getRV();
264 cluster_r[ncluster] = anode_mid_gap;
265 cluster_phi[ncluster] = phi;
266 cluster_phi_cc[ncluster] = phi_cc;
267 cluster_phi_tpc[ncluster] = phi_tpc;
268 cluster_ax_tpc[ncluster] = a_tpc;
269 cluster_bx_tpc[ncluster] = b_tpc;
272 cluster_r[ncluster] = anode_mid_gap;
273 cluster_v[ncluster] =
v;
274 cluster_v_cc[ncluster] = v_cc;
275 cluster_v_tpc[ncluster] = v_tpc;
276 cluster_av_tpc[ncluster] = a_tpc;
277 cluster_bv_tpc[ncluster] = b_tpc;
279 else if(flag==2 || flag==3) {
280 cluster_r[ncluster] = anode_mid_gap;
283 cluster_z[ncluster] = z;
284 cluster_z_cc[ncluster] = z_cc;
285 cluster_z_tpc[ncluster] = z_tpc;
286 cluster_phi[ncluster] = phi;
287 cluster_phi_cc[ncluster] = phi_cc;
288 cluster_phi_tpc[ncluster] = phi_tpc;
293 if(edep > tmp_charge_L1_S1) {
294 tmp_charge_L1_S1 = edep; tmp_cluster_L1_S1 = ncluster;
297 else if(layerid==1 && sheetid==0) {
298 if(edep > tmp_charge_L2_S1){
299 tmp_charge_L2_S1 = edep; tmp_cluster_L2_S1 = ncluster;
302 else if(layerid==1 && sheetid==1) {
303 if(edep > tmp_charge_L2_S2){
304 tmp_charge_L2_S2 = edep; tmp_cluster_L2_S2 = ncluster;
309 cluster_layerid[ncluster] = layerid;
310 cluster_sheetid[ncluster] = sheetid;
312 cluster_view[ncluster] = flag;
314 int flagB = (*iClusterCol)->getclusterflagb();
315 int flagE = (*iClusterCol)->getclusterflage();
317 if(flag==2 || flag==3) {
318 int clusterX = (*iClusterCol)->getclusterflagb();
319 RecCgemClusterCol::iterator iClusterCol2 = aClusterCol->begin()+clusterX;
320 int stripx1 = (*iClusterCol2)->getclusterflagb();
321 int stripx2 = (*iClusterCol2)->getclusterflage();
322 int sizex = stripx2 - stripx1 + 1;
323 double chargex = (*iClusterCol2)->getenergydeposit();
327 int clusterV = (*iClusterCol)->getclusterflage();
328 iClusterCol2 = aClusterCol->begin()+clusterV;
329 int stripv1 = (*iClusterCol2)->getclusterflagb();
330 int stripv2 = (*iClusterCol2)->getclusterflage();
331 int sizev = stripv2 - stripv1 + 1;
332 double chargev = (*iClusterCol2)->getenergydeposit();
341 cluster_idx[ncluster] = clusterX;
342 cluster_idv[ncluster] = clusterV;
343 cluster_stripx1[ncluster] = stripx1;
344 cluster_stripx2[ncluster] = stripx2;
345 cluster_stripv1[ncluster] = stripv1;
346 cluster_stripv2[ncluster] = stripv2;
347 cluster_sizex[ncluster] = sizex;
348 cluster_sizev[ncluster] = sizev;
349 cluster_qx[ncluster] = chargex;
350 cluster_qv[ncluster] = chargev;
351 cluster_ax_tpc[ncluster] = ax_tpc;
352 cluster_av_tpc[ncluster] = av_tpc;
353 cluster_bx_tpc[ncluster] = bx_tpc;
354 cluster_bv_tpc[ncluster] = bv_tpc;
357 int strip1 = (*iClusterCol)->getclusterflagb();
358 int strip2 = (*iClusterCol)->getclusterflage();
359 int size = strip2 - strip1 + 1;
360 double charge = (*iClusterCol)->getenergydeposit();
363 cluster_idx[ncluster] = clusterid;
364 cluster_sizex[ncluster] = size;
365 cluster_qx[ncluster] = charge;
366 cluster_stripx1[ncluster] = strip1;
367 cluster_stripx2[ncluster] = strip2;
370 cluster_idv[ncluster] = clusterid;
371 cluster_sizev[ncluster] = size;
372 cluster_qv[ncluster] = charge;
373 cluster_stripv1[ncluster] = strip1;
374 cluster_stripv2[ncluster] = strip2;
390 if(layerid == 0) ncluster_L1_S1_x++;
391 else if(layerid == 1) {
392 if(sheetid == 0) ncluster_L2_S1_x++;
393 else ncluster_L2_S2_x++;
398 if(layerid == 0) ncluster_L1_S1_v++;
399 else if(layerid == 1) {
400 if(sheetid == 0) ncluster_L2_S1_v++;
401 else ncluster_L2_S2_v++;
442 if(tmp_cluster_L1_S1!=-1) cluster_highest[tmp_cluster_L1_S1]=1;
443 if(tmp_cluster_L2_S1!=-1) cluster_highest[tmp_cluster_L2_S1]=1;
444 if(tmp_cluster_L2_S2!=-1) cluster_highest[tmp_cluster_L2_S2]=1;
451 return StatusCode::SUCCESS;
455 MsgStream log(
msgSvc(),name());
456 log << MSG::INFO <<
"TestCluster finalize()" << endreq;
459 return StatusCode::SUCCESS;
**********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
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
TestCluster(const std::string &name, ISvcLocator *pSvcLocator)