CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
TestClusterWithHit Class Reference

#include <TestClusterWithHit.h>

+ Inheritance diagram for TestClusterWithHit:

Public Member Functions

 TestClusterWithHit (const std::string &name, ISvcLocator *pSvcLocator)
 
 ~TestClusterWithHit ()
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 
void reset ()
 
void reset_hit ()
 
void reset_cluster_1d ()
 
void reset_cluster_2d ()
 
void reset_map ()
 

Detailed Description

Definition at line 24 of file TestClusterWithHit.h.

Constructor & Destructor Documentation

◆ TestClusterWithHit()

TestClusterWithHit::TestClusterWithHit ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Definition at line 44 of file TestClusterWithHit.cxx.

44 :
45 Algorithm(name,pSvcLocator){
46
47 // declareProperty("CosmicRayDataSetID", CosmicRayDataSetID = "CR201909");
48 declareProperty("selGoodDigi",m_selGoodDigi=1);
49 declareProperty("minDigiTime",m_minDigiTime=-8875);
50 declareProperty("maxDigiTime",m_maxDigiTime=-8562);
51 declareProperty("minQDigi",myQMin=0.0);
52 declareProperty("LUTfile", lutfile = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_from_17_to_17.root");
53
54
55
56}

◆ ~TestClusterWithHit()

TestClusterWithHit::~TestClusterWithHit ( )

Definition at line 58 of file TestClusterWithHit.cxx.

58 {
59 // delete m_SvcCgem;
60 // delete f;
61 // delete Tdigi;
62 // delete anode;
63 // delete tree;
64 // delete output;
65 delete lutreader;
66}

Member Function Documentation

◆ execute()

StatusCode TestClusterWithHit::execute ( )

read map std::map<int, std::vector < int > >::iterator it2; for(it2 = map_L1_S1_stripx_to_hit.begin(); it2 != map_L1_S1_stripx_to_hit.end(); it2++) { int istrip = it2->first; std::vector< int > hitlist = it2->second; cout << "L1/S1/x....reading strip " << istrip << " hit " << hitlist[0] << endl;
}

cout << "1D_CLUSTER index " << clusterid << " or " << ncluster << " layer " << layerid << " sheet " << sheetid << " view " << cluster_1d_view[ncluster] << endl; cout << " from " << flagB << " to " << flagE << endl; cout << "hits: ";

cout << "2D_CLUSTER index " << clusterid << " or " << ncluster << " layer " << layerid << " sheet " << sheetid << " view " << cluster_2d_view[ncluster] << endl; cout << " idx " << flagB << " idv " << flagE << endl;

cout << "************************" << endl; for(int iclu=0; iclu<ncluster; iclu++) { cout << "cluster 1d view " << cluster_1d_view[iclu] << " size " << cluster_1d_size[iclu] << " from " << cluster_1d_strip1[iclu] << " to " << cluster_1d_strip2[iclu] << endl;

for(int ihit = 0; ihit < cluster_1d_size[iclu] ; ihit++) { cout << "hit " << cluster_1d_hitindex[iclu][ihit] << endl; } }

Definition at line 303 of file TestClusterWithHit.cxx.

303 {
304
305 MsgStream log(msgSvc(), name());
306 if(event%1000==0) cout << "->TestClusterWithHit::execute " << event << endl;
307
308 //interface to event data service
309 ISvcLocator* svcLocator = Gaudi::svcLocator();
310 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
311 if (sc.isFailure())
312 cout<<"Could not accesss EventDataSvc!"<<endl;
313
314 // reset
315 reset();
316 // -------------------------
317
318 //retrieve CGEM hits from TDS
319 SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
320 if(!aDigiCol)
321 cout<<"Could not retrieve CGEM digi collection"<<endl;
322 nhit = aDigiCol->size();
323
324 // loop on hits
325 int ihit=0;
326
327 CgemDigiCol::iterator iDigiCol;
328 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
329 {
330 const Identifier ident = (*iDigiCol)->identify();
331 int strip = CgemID::strip(ident);
332 int layer = CgemID::layer(ident);
333 int sheet = CgemID::sheet(ident);
334 int view = 1;
335 bool is_xstrip = CgemID::is_xstrip(ident);
336 if(is_xstrip == true) view = 0;
337 double charge = (*iDigiCol)->getCharge_fc();
338 double time = (*iDigiCol)->getTime_ns();
339
340 hit_strip[ihit] = strip;
341 hit_layer[ihit] = layer;
342 hit_sheet[ihit] = sheet;
343 hit_view[ihit] = view;
344
345 if(selDigi(iDigiCol, m_selGoodDigi)==true) {
346
347 // cout << "HIT " << ihit << " strip " << strip << " layer " << layer << " sheet " << sheet << " view " << view << endl;
348 // map strips to hits (maybe more hits for the same strip)
349 std::map<int, std::vector < int > > *map_strip_to_hit = GetStripToHitMap(layer, sheet, view);
350 std::map<int, std::vector < int > >::iterator it1 = map_strip_to_hit->find(strip);
351 if (it1 != map_strip_to_hit->end()) {
352 it1->second.push_back(ihit);
353 }
354 else {
355 std::vector< int > hitlist;
356 hitlist.push_back(ihit);
357 std::pair<int, std::vector < int> > pair(strip, hitlist);
358 map_strip_to_hit->insert(pair);
359 }
360 }
361 // --------------------------------------------------------
362
363 ihit++;
364 }
365 /**
366 //read map
367 std::map<int, std::vector < int > >::iterator it2;
368 for(it2 = map_L1_S1_stripx_to_hit.begin(); it2 != map_L1_S1_stripx_to_hit.end(); it2++)
369 {
370 int istrip = it2->first;
371 std::vector< int > hitlist = it2->second;
372 cout << "L1/S1/x....reading strip " << istrip << " hit " << hitlist[0] << endl;
373 }
374 **/
375
376 // -------------------------------------------------------------------
377 //retrieve CGEM clusters from TDS
378 SmartDataPtr<RecCgemClusterCol> aClusterCol(m_evtSvc,"/Event/Recon/RecCgemClusterCol");
379 if(!aClusterCol)
380 cout<<"Could not retrieve CGEM cluster collection"<<endl;
381
382 int nclu = aClusterCol->size();
383 if(nclu > MAXNOFCLUSTERS) {
384 std::cout << "nclu " << nclu << std::endl;
385 nclu = MAXNOFCLUSTERS; // CHECK HACK
386 }
387
388 int nclusterx = 0;
389 int nclusterv = 0;
390 int nclusterxv = 0;
391
392 anode = m_SvcCgem->getReadoutPlane(0, 0);
393 anode_radius_L1_x = anode->getRX();
394 anode_radius_L1_v = anode->getRV();
395 anode_mid_gap_L1 = anode->getMidRAtGap();
396 anode = m_SvcCgem->getReadoutPlane(1, 0);
397 anode_radius_L2_x = anode->getRX();
398 anode_radius_L2_v = anode->getRV();
399 anode_mid_gap_L2 = anode->getMidRAtGap();
400
401 int tmp_cluster_L1_S1 = -1;
402 int tmp_cluster_L2_S1 = -1;
403 int tmp_cluster_L2_S2 = -1;
404 double tmp_charge_L1_S1 = 0.;
405 double tmp_charge_L2_S1 = 0.;
406 double tmp_charge_L2_S2 = 0.;
407
408 ncluster = 0;
409 RecCgemClusterCol::iterator iClusterCol;
410 for(iClusterCol=aClusterCol->begin(); iClusterCol!=aClusterCol->end(); iClusterCol++)
411 {
412
413 if(ncluster==MAXNOFCLUSTERS) break;
414
415 int flag = (*iClusterCol)->getflag(); // 0=x 1=v 2=xv
416 int clusterid = (*iClusterCol)->getclusterid();
417 int trkid = (*iClusterCol)->getTrkId();
418 int layerid = (*iClusterCol)->getlayerid();
419 int sheetid = (*iClusterCol)->getsheetid();
420 double edep = (*iClusterCol)->getenergydeposit();
421 double phi = (*iClusterCol)->getrecphi();
422 double v = (*iClusterCol)->getrecv();
423 double z = (*iClusterCol)->getRecZ();
424 double phi_cc = (*iClusterCol)->getrecphi_CC();
425 double v_cc = (*iClusterCol)->getrecv_CC();
426 double z_cc = (*iClusterCol)->getRecZ_CC();
427 double phi_tpc = (*iClusterCol)->getrecphi_mTPC();
428 double v_tpc = (*iClusterCol)->getrecv_mTPC();
429 double z_tpc = (*iClusterCol)->getRecZ_mTPC();
430 double a_tpc = 0; //(*iClusterCol)->getparamA_mTPC();
431 double b_tpc = 0; //(*iClusterCol)->getparamB_mTPC();
432
433 if(sheetid == -1) {
434 ncluster++;
435 continue; // CHECK HACK FOR NOW
436 }
437
438 // anode mid gap radius
439 anode = m_SvcCgem->getReadoutPlane(layerid, sheetid);
440 double anode_mid_gap = anode->getMidRAtGap();
441
442 // beginning and ending flags
443 // if 2D cluster --> they are the x & v 1D clusters
444 // if 1D cluster --> they are the first and last contigous strips
445 int flagB = (*iClusterCol)->getclusterflagb();
446 int flagE = (*iClusterCol)->getclusterflage();
447
448 // ------------------------------------------------------------------- 1D cluster
449 // positions CC & uTPC
450 if(flag==0 || flag == 1) {
451 cluster_1d_ID[ncluster] = clusterid;
452 // time & charge
453 cluster_1d_t[ncluster] = 0; // CHECK
454 cluster_1d_q[ncluster] = edep;
455 // layer/sheet/view
456 cluster_1d_layerid[ncluster] = layerid;
457 cluster_1d_sheetid[ncluster] = sheetid;
458 cluster_1d_view[ncluster] = flag;
459 cluster_1d_r[ncluster] = anode_mid_gap;
460 cluster_1d_phi[ncluster] = phi;
461 cluster_1d_phi_cc[ncluster] = phi_cc;
462 cluster_1d_phi_tpc[ncluster] = phi_tpc;
463 cluster_1d_v[ncluster] = v;
464 cluster_1d_v_cc[ncluster] = v_cc;
465 cluster_1d_v_tpc[ncluster] = v_tpc;
466 cluster_1d_a_tpc[ncluster] = a_tpc;
467 cluster_1d_b_tpc[ncluster] = b_tpc;
468 // from strip1 to strip2
469 cluster_1d_strip1[ncluster] = flagB;
470 cluster_1d_strip2[ncluster] = flagE;
471 int size = flagE - flagB + 1;
472 cluster_1d_size[ncluster] = size;
473
474 // cout << "getting the hits from strip1 "<< flagB << " to strip2 " << flagE << " size " << size << endl;
475 // get the hits
476 std::map<int, std::vector < int > > *read_map_strip_to_hit = GetStripToHitMap(layerid, sheetid, flag);
477 std::map<int, std::vector < int > >::iterator it2;
478 /**
479 cout << "1D_CLUSTER index " << clusterid << " or " << ncluster
480 << " layer " << layerid << " sheet " << sheetid
481 << " view " << cluster_1d_view[ncluster] << endl;
482 cout << " from " << flagB << " to " << flagE << endl;
483 cout << "hits: ";
484 **/
485
486 // if(size >= 50) cout << "size " << size << endl;
487 for(int istrip = 0; istrip < size; istrip++) {
488 int stripid = flagB + istrip;
489 it2 = read_map_strip_to_hit->find(stripid);
490 std::vector< int > hitlist = it2->second;
491 int size_hitlist = hitlist.size();
492 // if(size_hitlist>1) cout << event << " cluster " << ncluster << " hitlist > 1 " << size_hitlist << endl;
493 cluster_1d_hitindex[ncluster][istrip] = hitlist[size_hitlist-1]; // take only the first hit CHECK THIS HAS TO BE CHANGED!
494 // cout << " " << hitlist[0];
495 }
496 // cout << endl;
497 // -----------
498 ncluster_1d++;
499 }
500 else if(flag==2 || flag==3) { // ------------------------------------- 2D - xv view
501 cluster_2d_ID[ncluster] = clusterid;
502 // time & charge
503 cluster_2d_t[ncluster] = 0; // CHECK
504 cluster_2d_q[ncluster] = edep;
505 // layer/sheet/view
506 cluster_2d_layerid[ncluster] = layerid;
507 cluster_2d_sheetid[ncluster] = sheetid;
508 cluster_2d_view[ncluster] = flag;
509 cluster_2d_r[ncluster] = anode_mid_gap;
510 cluster_2d_z[ncluster] = z;
511 cluster_2d_z_cc[ncluster] = z_cc;
512 cluster_2d_z_tpc[ncluster] = z_tpc;
513 cluster_2d_phi[ncluster] = phi;
514 cluster_2d_phi_cc[ncluster] = phi_cc;
515 cluster_2d_phi_tpc[ncluster] = phi_tpc;
516 // composing 1D_clusters
517 cluster_2d_idx[ncluster] = flagB;
518 cluster_2d_idv[ncluster] = flagE;
519
520 /** cout << "2D_CLUSTER index " << clusterid << " or " << ncluster
521 << " layer " << layerid << " sheet " << sheetid
522 << " view " << cluster_2d_view[ncluster] << endl;
523 cout << " idx " << flagB << " idv " << flagE << endl;
524
525 **/
526 // find the highest charge cluster
527 if(flag==2) { // CHECK THIS
528 if(layerid==0) {
529 if(edep > tmp_charge_L1_S1) {
530 tmp_charge_L1_S1 = edep; tmp_cluster_L1_S1 = ncluster;
531 }
532 }
533 else if(layerid==1 && sheetid==0) {
534 if(edep > tmp_charge_L2_S1){
535 tmp_charge_L2_S1 = edep; tmp_cluster_L2_S1 = ncluster;
536 }
537 }
538 else if(layerid==1 && sheetid==1) {
539 if(edep > tmp_charge_L2_S2){
540 tmp_charge_L2_S2 = edep; tmp_cluster_L2_S2 = ncluster;
541 }
542 }
543 }
544 // -----------
545 ncluster_2d++;
546 }
547
548 // cluster counters
549 ncluster++;
550 if(flag == 0) {
551 nclusterx++;
552 if(layerid == 0) ncluster_1d_L1_S1_x++;
553 else if(layerid == 1) {
554 if(sheetid == 0) ncluster_1d_L2_S1_x++;
555 else ncluster_1d_L2_S2_x++;
556 }
557 }
558 else if(flag == 1) {
559 nclusterv++;
560 if(layerid == 0) ncluster_1d_L1_S1_v++;
561 else if(layerid == 1) {
562 if(sheetid == 0) ncluster_1d_L2_S1_v++;
563 else ncluster_1d_L2_S2_v++;
564 }
565 }
566 else if(flag == 2) {
567 nclusterxv++;
568 if(layerid == 0) ncluster_2d_L1_S1++;
569 else if(layerid == 1) {
570 if(sheetid == 0) ncluster_2d_L2_S1++;
571 else ncluster_2d_L2_S2++;
572 }
573 }
574 // --------
575 }
576
577 // std::cout << "cluster " << tmp_cluster_L1_S1 << " " << tmp_cluster_L2_S1 << " " << tmp_cluster_L2_S2 << std::endl;
578 // std::cout << "charge " << tmp_charge_L1_S1 << " " << tmp_charge_L2_S1 << " " << tmp_charge_L2_S2 << std::endl;
579
580 // set the highest charge 2D - cluster
581 if(tmp_cluster_L1_S1!=-1) cluster_2d_highest[tmp_cluster_L1_S1]=1;
582 if(tmp_cluster_L2_S1!=-1) cluster_2d_highest[tmp_cluster_L2_S1]=1;
583 if(tmp_cluster_L2_S2!=-1) cluster_2d_highest[tmp_cluster_L2_S2]=1;
584
585 /**
586 cout << "************************" << endl;
587 for(int iclu=0; iclu<ncluster; iclu++) {
588 cout << "cluster 1d view " << cluster_1d_view[iclu] << " size " << cluster_1d_size[iclu] << " from " << cluster_1d_strip1[iclu] << " to " << cluster_1d_strip2[iclu] << endl;
589
590 for(int ihit = 0; ihit < cluster_1d_size[iclu] ; ihit++) {
591 cout << "hit " << cluster_1d_hitindex[iclu][ihit] << endl;
592 }
593 }
594 **/
595
596 tree->Fill();
597 event++;
598
599 return StatusCode::SUCCESS;
600}
time
Double_t time
Definition: DataBase/tau_mode.c:7
v
**********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
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
MAXNOFCLUSTERS
#define MAXNOFCLUSTERS
Definition: TestClusterWithHit.h:21
CgemID::strip
static int strip(const Identifier &id)
Definition: CgemID.cxx:83
CgemID::sheet
static int sheet(const Identifier &id)
Definition: CgemID.cxx:77
CgemID::layer
static int layer(const Identifier &id)
Definition: CgemID.cxx:71
CgemID::is_xstrip
static bool is_xstrip(const Identifier &id)
Definition: CgemID.cxx:64
ICgemGeomSvc::getReadoutPlane
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0

◆ finalize()

StatusCode TestClusterWithHit::finalize ( )

Definition at line 602 of file TestClusterWithHit.cxx.

602 {
603 MsgStream log(msgSvc(),name());
604 log << MSG::INFO << "TestClusterWithHit finalize()" << endreq;
605 output->Write();
606 output->Close();
607 return StatusCode::SUCCESS;
608}

◆ initialize()

StatusCode TestClusterWithHit::initialize ( )

Definition at line 68 of file TestClusterWithHit.cxx.

68 {
69 MsgStream log(msgSvc(), name());
70 log << MSG::INFO << "TestClusterWithHit initialize()" << endreq;
71
72
73 // CgemGeomSvc
74 StatusCode sc;
75 sc = service("CgemGeomSvc", m_SvcCgem);
76 if(sc != StatusCode::SUCCESS) {
77 log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
78 return StatusCode::FAILURE;
79 }
80
81 // LUT
82 lutreader = new CgemLUTReader(lutfile);
83 lutreader->ReadLUT();
84
85 output = new TFile("cluster.root", "RECREATE");
86 tree = new TTree("tree","cluster info tree");
87 DefineClusterTree();
88 DefineHitTree();
89
90 event=0;
91
92 return StatusCode::SUCCESS;
93
94}

◆ reset()

void TestClusterWithHit::reset ( )

Definition at line 159 of file TestClusterWithHit.cxx.

159 {
160 reset_map();
161 reset_cluster_1d();
162 reset_cluster_2d();
163 reset_hit();
164
165}

Referenced by execute().

◆ reset_cluster_1d()

void TestClusterWithHit::reset_cluster_1d ( )

Definition at line 178 of file TestClusterWithHit.cxx.

178 {
179
180 ncluster_1d = 0;
181 ncluster_1d_L1_S1_x = 0;
182 ncluster_1d_L2_S1_x = 0;
183 ncluster_1d_L2_S2_x = 0;
184 ncluster_1d_L1_S1_v = 0;
185 ncluster_1d_L2_S1_v = 0;
186 ncluster_1d_L2_S2_v = 0;
187
188 for(int iclu = 0; iclu < MAXNOFCLUSTERS; iclu++)
189 {
190 cluster_1d_t[iclu] = 0.;
191 cluster_1d_q[iclu] = 0.;
192 cluster_1d_r[iclu] = 0.;
193 cluster_1d_v[iclu] = 0.;
194 cluster_1d_phi[iclu] = 0.;
195 cluster_1d_v_cc[iclu] = 0.;
196 cluster_1d_phi_cc[iclu] = 0.;
197 cluster_1d_v_tpc[iclu] = 0.;
198 cluster_1d_phi_tpc[iclu] = 0.;
199 cluster_1d_a_tpc[iclu] = 0.;
200 cluster_1d_b_tpc[iclu] = 0.;
201 cluster_1d_layerid[iclu] = -1;
202 cluster_1d_sheetid[iclu] = -1;
203 cluster_1d_view[iclu] = -1;
204 cluster_1d_size[iclu] = -1;
205 cluster_1d_strip1[iclu] = -1;
206 cluster_1d_strip2[iclu] = -1;
207 for(int ihit = 0; ihit < MAXNOFHITS; ihit++) cluster_1d_hitindex[iclu][ihit] = -1;
208 }
209
210
211
212}
MAXNOFHITS
#define MAXNOFHITS
Definition: TestClusterWithHit.h:22

Referenced by reset().

◆ reset_cluster_2d()

void TestClusterWithHit::reset_cluster_2d ( )

Definition at line 214 of file TestClusterWithHit.cxx.

214 {
215
216 ncluster_2d = 0;
217 ncluster_2d_L1_S1 = 0;
218 ncluster_2d_L2_S1 = 0;
219 ncluster_2d_L2_S2 = 0;
220
221 for(int iclu = 0; iclu < MAXNOFCLUSTERS; iclu++)
222 {
223 cluster_2d_t[iclu] = 0.;
224 cluster_2d_q[iclu] = 0.;
225 cluster_2d_r[iclu] = 0.;
226 cluster_2d_z[iclu] = 0.;
227 cluster_2d_phi[iclu] = 0.;
228 cluster_2d_z_cc[iclu] = 0.;
229 cluster_2d_phi_cc[iclu] = 0.;
230 cluster_2d_z_tpc[iclu] = 0.;
231 cluster_2d_phi_tpc[iclu] = 0.;
232 cluster_2d_layerid[iclu] = -1;
233 cluster_2d_sheetid[iclu] = -1;
234 cluster_2d_view[iclu] = -1;
235 cluster_2d_idx[iclu] = -1;
236 cluster_2d_idv[iclu] = -1;
237 cluster_2d_highest[iclu] = -1;
238 }
239}

Referenced by reset().

◆ reset_hit()

void TestClusterWithHit::reset_hit ( )

Definition at line 166 of file TestClusterWithHit.cxx.

166 {
167 // event = 0;
168 nhit = 0;
169
170 for(int ihit = 0; ihit < MAXNOFHITS; ihit++) {
171 hit_strip[ihit] = -1;
172 hit_view[ihit] = -1;
173 hit_layer[ihit] = -1;
174 hit_sheet[ihit] = -1;
175 }
176}

Referenced by reset().

◆ reset_map()

void TestClusterWithHit::reset_map ( )

Definition at line 243 of file TestClusterWithHit.cxx.

243 {
244 map_L1_S1_stripx_to_hit.clear();
245 map_L2_S1_stripx_to_hit.clear();
246 map_L2_S2_stripx_to_hit.clear();
247 map_L1_S1_stripv_to_hit.clear();
248 map_L2_S1_stripv_to_hit.clear();
249 map_L2_S2_stripv_to_hit.clear();
250}

Referenced by reset().

The documentation for this class was generated from the following files: