CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
TestCluster.cxx
Go to the documentation of this file.
1// **************************************************************************/
2// authors: L. Lavezzi (univ. of Torino & INFN, Italy)
3//
4
5//include system lib
6#include <iostream>
7#include <iomanip>
8#include <string>
9#include <cmath>
10
11// Include files
12#include "GaudiKernel/AlgFactory.h"
13#include "GaudiKernel/DataObject.h"
14#include "GaudiKernel/IEventProcessor.h"
15
16#include "GaudiKernel/Incident.h"
17#include "GaudiKernel/IIncidentSvc.h"
18#include "GaudiKernel/Memory.h"
19
20#include <csignal>
21
22//for save digit & cluster
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"
28
31#include "Identifier/CgemID.h"
32
34#include "RawEvent/DigiEvent.h"
37#include "GaudiKernel/SmartDataPtr.h"
38
40
41#include "TMath.h"
42//using namespace std;
43
44TestCluster::TestCluster(const std::string& name, ISvcLocator* pSvcLocator):
45 Algorithm(name,pSvcLocator){
46
47 // declareProperty("CosmicRayDataSetID", CosmicRayDataSetID = "CR201909");
48
49
50}
51
53 // delete m_SvcCgem;
54 // delete f;
55 // delete Tdigi;
56 // delete anode;
57 // delete tree;
58 // delete output;
59
60}
61
63 MsgStream log(msgSvc(), name());
64 log << MSG::INFO << "TestCluster initialize()" << endreq;
65
66
67 // CgemGeomSvc
68 StatusCode sc;
69 sc = service("CgemGeomSvc", m_SvcCgem);
70 if(sc != StatusCode::SUCCESS) {
71 log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
72 return StatusCode::FAILURE;
73 }
74
75 output = new TFile("cluster.root", "RECREATE");
76 tree = new TTree("tree","cluster info tree");
77
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");
121
122 return StatusCode::SUCCESS;
123}
124
126 ncluster = 0;
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;
133
134 for(int iclu = 0; iclu < MAXNOFCLUSTERS; iclu++)
135 {
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;
166 }
167}
168
170
171 MsgStream log(msgSvc(), name());
172 // cout << "->TestCluster::execute" << endl;
173
174 //interface to event data service
175 ISvcLocator* svcLocator = Gaudi::svcLocator();
176 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
177 if (sc.isFailure())
178 cout<<"Could not accesss EventDataSvc!"<<endl;
179
180 //retrieve CGEM clusters from TDS
181 SmartDataPtr<RecCgemClusterCol> aClusterCol(m_evtSvc,"/Event/Recon/RecCgemClusterCol");
182 if(!aClusterCol)
183 cout<<"Could not retrieve CGEM cluster collection"<<endl;
184
185 int nclu = aClusterCol->size();
186
187 if(nclu > MAXNOFCLUSTERS) {
188 std::cout << "nclu " << nclu << std::endl;
189 nclu = MAXNOFCLUSTERS; // CHECK HACK
190 }
191
192 reset();
193
194 int nclusterx = 0;
195 int nclusterv = 0;
196 int nclusterxv = 0;
197
198 anode = m_SvcCgem->getReadoutPlane(0, 0);
199 anode_radius_L1_x = anode->getRX();
200 anode_radius_L1_v = anode->getRV();
201 anode_mid_gap_L1 = anode->getMidRAtGap();
202 anode = m_SvcCgem->getReadoutPlane(1, 0);
203 anode_radius_L2_x = anode->getRX();
204 anode_radius_L2_v = anode->getRV();
205 anode_mid_gap_L2 = anode->getMidRAtGap();
206
207
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.;
214
215 RecCgemClusterCol::iterator iClusterCol;
216 for(iClusterCol=aClusterCol->begin(); iClusterCol!=aClusterCol->end(); iClusterCol++)
217 {
218
219 if(ncluster==MAXNOFCLUSTERS) break;
220
221 int flag = (*iClusterCol)->getflag(); // 0=x 1=v 2=xv
222
223 // cout << flag << " " << (*iClusterCol)->getrecphi() << " " << (*iClusterCol)->getrecv() << " " << (*iClusterCol)->getRecZ() << endl;
224 if(flag==2 || flag==3) nclusterxv++;
225
226 int clusterid = (*iClusterCol)->getclusterid();
227 int trkid = (*iClusterCol)->getTrkId();
228 int layerid = (*iClusterCol)->getlayerid();
229 int sheetid = (*iClusterCol)->getsheetid();
230
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();
241 double a_tpc = 0; //(*iClusterCol)->getparamA_mTPC();
242 double b_tpc = 0; //(*iClusterCol)->getparamB_mTPC();
243
244 // if(edep==0) continue; // CHECK HACK FOR NOW
245 if(sheetid==-1) {
246 // cout << "TestCluster: SHEETID " << sheetid << " on layer " << layerid << " and view " << flag << endl;
247 continue; // CHECK HACK FOR NOW
248 }
249 // cout << ncluster << " A " << flag << " " << phi << " " << v << endl;
250
251 cluster_t[ncluster] = 0; // CHECK
252 // cout << "edep " << edep << endl;
253 cluster_q[ncluster] = edep;
254 // cout << "layerid " <<layerid << endl;
255 // cout << "sheetid " << sheetid << endl;
256 anode = m_SvcCgem->getReadoutPlane(layerid, sheetid);
257 double anode_radius_x = anode->getRX();
258 double anode_radius_v = anode->getRV();
259 double anode_mid_gap = anode->getMidRAtGap();
260 // cout << ncluster << " B " << " flag " << flag << " phi " << phi << " v " << v << " z " << z << endl;
261 if(flag==0) {
262 // cluster_x[ncluster] = anode_radius_x * TMath::Cos(phi);
263 // cluster_y[ncluster] = anode_radius_x * TMath::Sin(phi);
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;
270 }
271 else if(flag==1) {
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;
278 }
279 else if(flag==2 || flag==3) {
280 cluster_r[ncluster] = anode_mid_gap;
281 // cluster_x[ncluster] = 0.5* (anode_radius_x + anode_radius_v) * TMath::Cos(phi);
282 // cluster_y[ncluster] = 0.5* (anode_radius_x + anode_radius_v) * TMath::Sin(phi);
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;
289 // cout << "flag " << flag << " layerid " << layerid << " sheetid " << sheetid << " charge " << edep
290 // << " ncluster " << ncluster << " indice " << iClusterCol-aClusterCol->begin() << endl;
291
292 if(layerid==0) {
293 if(edep > tmp_charge_L1_S1) {
294 tmp_charge_L1_S1 = edep; tmp_cluster_L1_S1 = ncluster;
295 }
296 }
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;
300 }
301 }
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;
305 }
306 }
307 }
308
309 cluster_layerid[ncluster] = layerid;
310 cluster_sheetid[ncluster] = sheetid;
311
312 cluster_view[ncluster] = flag;
313
314 int flagB = (*iClusterCol)->getclusterflagb();
315 int flagE = (*iClusterCol)->getclusterflage();
316
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();
324 double ax_tpc = 0; // (*iClusterCol2)->getparamA_mTPC();
325 double bx_tpc = 0; // ((*iClusterCol2)->getparamB_mTPC();
326
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();
333 double av_tpc = 0; // ((*iClusterCol2)->getparamA_mTPC();
334 double bv_tpc = 0; // ((*iClusterCol2)->getparamB_mTPC();
335
336 // cout << "cluster index " << iClusterCol-aClusterCol->begin()
337 // << " view " << cluster_view[ncluster]
338 // << " cluster x " << clusterX << " from " << stripx1 << " to " << stripx2
339 // << " cluster v " << clusterV << " from " << stripv1 << " to " << stripv2 << endl;
340
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;
355 }
356 else {
357 int strip1 = (*iClusterCol)->getclusterflagb();
358 int strip2 = (*iClusterCol)->getclusterflage();
359 int size = strip2 - strip1 + 1;
360 double charge = (*iClusterCol)->getenergydeposit();
361
362 if(flag==0) {
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;
368 }
369 else if(flag==1) {
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;
375 }
376
377
378 // cout << "cluster index " << iClusterCol-aClusterCol->begin()
379 // << " view " << cluster_view[ncluster]
380 // << " strip1 " << strip1
381 // << " strip2 " << strip2 << endl;
382
383 }
384
385
386
387 ncluster++;
388 if(flag == 0) {
389 nclusterx++;
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++;
394 }
395 }
396 else if(flag == 1) {
397 nclusterv++;
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++;
402 }
403 }
404
405
406 // int getclusterflagb(void) const{return m_clusterflag.first;};
407 // int getclusterflage(void) const{return m_clusterflag.second;};
408 // std::vector<int> getStripsMTPC() {return m_strips_mTPC;};
409
410
411 // if(flag==0) {
412 // cout<< "X VIEW" << endl;
413 // cout << "clusterid " << clusterid << " edep " << edep << endl;
414 // cout << "layerid " << layerid << " sheetid " << sheetid << endl;
415 // cout << "phi " << phi << " cc " << phi_cc << " tpc " << phi_tpc << endl;
416 // cout << "v " << v << " cc " << v_cc << " tpc " << v_tpc << endl;
417 // cout << "z " << z << " cc " << z_cc << " tpc " << z_tpc << endl;
418 // }
419 // else if(flag==1) {
420 // cout<< "V VIEW" << endl;
421 // cout << "clusterid " << clusterid << " edep " << edep << endl;
422 // cout << "layerid " << layerid << " sheetid " << sheetid << endl;
423 // cout << "phi " << phi << " cc " << phi_cc << " tpc " << phi_tpc << endl;
424 // cout << "v " << v << " cc " << v_cc << " tpc " << v_tpc << endl;
425 // cout << "z " << z << " cc " << z_cc << " tpc " << z_tpc << endl;
426 // }
427 // else if(flag==2) {
428 // cout<< "XV VIEW" << endl;
429 // cout << "clusterid " << clusterid << " edep " << edep << endl;
430 // cout << "layerid " << layerid << " sheetid " << sheetid << endl;
431 // cout << "phi " << phi << " cc " << phi_cc << " tpc " << phi_tpc << endl;
432 // cout << "v " << v << " cc " << v_cc << " tpc " << v_tpc << endl;
433 // cout << "z " << z << " cc " << z_cc << " tpc " << z_tpc << endl;
434 // }
435
436
437 }
438
439 // std::cout << "cluster " << tmp_cluster_L1_S1 << " " << tmp_cluster_L2_S1 << " " << tmp_cluster_L2_S2 << std::endl;
440 // std::cout << "charge " << tmp_charge_L1_S1 << " " << tmp_charge_L2_S1 << " " << tmp_charge_L2_S2 << std::endl;
441
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;
445
446 tree->Fill();
447
448
449
450
451 return StatusCode::SUCCESS;
452}
453
455 MsgStream log(msgSvc(),name());
456 log << MSG::INFO << "TestCluster finalize()" << endreq;
457 output->Write();
458 output->Close();
459 return StatusCode::SUCCESS;
460}
461
462
**********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
Definition KarLud.h:35
IMessageSvc * msgSvc()
#define MAXNOFCLUSTERS
Definition TestCluster.h:18
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
StatusCode initialize()
StatusCode finalize()
TestCluster(const std::string &name, ISvcLocator *pSvcLocator)
StatusCode execute()