CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
TestClusterWithHit.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
29#include "CgemRawEvent/CgemDigi.h"
30#include "CgemRecEvent/RecCgemCluster.h"
31#include "Identifier/CgemID.h"
32
33#include "RawEvent/RawDataUtil.h"
34#include "RawEvent/DigiEvent.h"
35#include "ReconEvent/ReconEvent.h"
36#include "EventModel/EventHeader.h"
37#include "GaudiKernel/SmartDataPtr.h"
38
39#include "CgemClusterCreate/TestClusterWithHit.h"
40
41#include "TMath.h"
42//using namespace std;
43
44TestClusterWithHit::TestClusterWithHit(const std::string& name, ISvcLocator* pSvcLocator):
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}
57
58TestClusterWithHit::~TestClusterWithHit(){
59 // delete m_SvcCgem;
60 // delete f;
61 // delete Tdigi;
62 // delete anode;
63 // delete tree;
64 // delete output;
65 delete lutreader;
66}
67
68StatusCode TestClusterWithHit::initialize(){
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}
95
96void TestClusterWithHit::DefineHitTree() {
97
98 tree->Branch("event", &event, "event/I");
99 tree->Branch("nhit", &nhit, "nhit/I");
100 // from geometry
101 tree->Branch("hit_strip", &hit_strip, "hit_strip[nhit]/I");
102 tree->Branch("hit_view", &hit_view, "hit_view[nhit]/I");
103 tree->Branch("hit_layer", &hit_layer, "hit_layer[nhit]/I");
104 tree->Branch("hit_sheet", &hit_sheet, "hit_sheet[nhit]/I");
105
106
107}
108void TestClusterWithHit::DefineClusterTree() {
109
110 tree->Branch("ncluster", &ncluster,"ncluster/I");
111 // 1d
112 tree->Branch("ncluster_1d", &ncluster_1d,"ncluster_1d/I");
113 tree->Branch("ncluster_1d_L1_S1_x", &ncluster_1d_L1_S1_x,"ncluster_1d_L1_S1_x/I");
114 tree->Branch("ncluster_1d_L2_S1_x", &ncluster_1d_L2_S1_x,"ncluster_1d_L2_S1_x/I");
115 tree->Branch("ncluster_1d_L2_S2_x", &ncluster_1d_L2_S2_x,"ncluster_1d_L2_S2_x/I");
116 tree->Branch("ncluster_1d_L1_S1_v", &ncluster_1d_L1_S1_v,"ncluster_1d_L1_S1_v/I");
117 tree->Branch("ncluster_1d_L2_S1_v", &ncluster_1d_L2_S1_v,"ncluster_1d_L2_S1_v/I");
118 tree->Branch("ncluster_1d_L2_S2_v", &ncluster_1d_L2_S2_v,"ncluster_1d_L2_S2_v/I");
119 tree->Branch("cluster_1d_ID", &cluster_1d_ID,"cluster_1d_ID[ncluster]/D");
120 tree->Branch("cluster_1d_t", &cluster_1d_t,"cluster_1d_t[ncluster]/D");
121 tree->Branch("cluster_1d_q", &cluster_1d_q,"cluster_1d_q[ncluster]/D");
122 tree->Branch("cluster_1d_r", &cluster_1d_r,"cluster_1d_r[ncluster]/D");
123 tree->Branch("cluster_1d_v", &cluster_1d_v, "cluster_1d_v[ncluster]/D");
124 tree->Branch("cluster_1d_v_cc", &cluster_1d_v_cc, "cluster_1d_v_cc[ncluster]/D");
125 tree->Branch("cluster_1d_v_tpc", &cluster_1d_v_tpc, "cluster_1d_v_tpc[ncluster]/D");
126 tree->Branch("cluster_1d_phi", &cluster_1d_phi, "cluster_1d_phi[ncluster]/D");
127 tree->Branch("cluster_1d_phi_cc", &cluster_1d_phi_cc, "cluster_1d_phi_cc[ncluster]/D");
128 tree->Branch("cluster_1d_phi_tpc", &cluster_1d_phi_tpc, "cluster_1d_phi_tpc[ncluster]/D");
129 tree->Branch("cluster_1d_layerid", &cluster_1d_layerid,"cluster_1d_layerid[ncluster]/I");
130 tree->Branch("cluster_1d_sheetid", &cluster_1d_sheetid,"cluster_1d_sheetid[ncluster]/I");
131 tree->Branch("cluster_1d_view", &cluster_1d_view,"cluster_1d_view[ncluster]/I");
132 tree->Branch("cluster_1d_strip1", &cluster_1d_strip1, "cluster_1d_strip1[ncluster]/I");
133 tree->Branch("cluster_1d_strip2", &cluster_1d_strip2, "cluster_1d_strip2[ncluster]/I");
134 tree->Branch("cluster_1d_size", &cluster_1d_size,"cluster_1d_size[ncluster]/I");
135 tree->Branch("cluster_1d_a_tpc", &cluster_1d_a_tpc,"cluster_1d_a_tpc[ncluster]/D");
136 tree->Branch("cluster_1d_b_tpc", &cluster_1d_b_tpc,"cluster_1d_b_tpc[ncluster]/D");
137 tree->Branch("cluster_1d_hitindex", cluster_1d_hitindex,"cluster_1d_hitindex[ncluster][50]/I");
138
139 // 2d
140 tree->Branch("ncluster_2d", &ncluster_2d,"ncluster_2d/I");
141 tree->Branch("cluster_2d_t", &cluster_2d_t,"cluster_2d_t[ncluster]/D");
142 tree->Branch("cluster_2d_q", &cluster_2d_q,"cluster_2d_q[ncluster]/D");
143 tree->Branch("cluster_2d_r", &cluster_2d_r,"cluster_2d_r[ncluster]/D");
144 tree->Branch("cluster_2d_z", &cluster_2d_z,"cluster_2d_z[ncluster]/D");
145 tree->Branch("cluster_2d_z_cc", &cluster_2d_z_cc, "cluster_2d_z_cc[ncluster]/D");
146 tree->Branch("cluster_2d_z_tpc", &cluster_2d_z_tpc, "cluster_2d_z_tpc[ncluster]/D");
147 tree->Branch("cluster_2d_phi", &cluster_2d_phi, "cluster_2d_phi[ncluster]/D");
148 tree->Branch("cluster_2d_phi_cc", &cluster_2d_phi_cc, "cluster_2d_phi_cc[ncluster]/D");
149 tree->Branch("cluster_2d_phi_tpc", &cluster_2d_phi_tpc, "cluster_2d_phi_tpc[ncluster]/D");
150 tree->Branch("cluster_2d_layerid", &cluster_2d_layerid,"cluster_2d_layerid[ncluster]/I");
151 tree->Branch("cluster_2d_sheetid", &cluster_2d_sheetid,"cluster_2d_sheetid[ncluster]/I");
152 tree->Branch("cluster_2d_view", &cluster_2d_view,"cluster_2d_view[ncluster]/I");
153 tree->Branch("cluster_2d_idx", &cluster_2d_idx, "cluster_2d_idx[ncluster]/I");
154 tree->Branch("cluster_2d_idv", &cluster_2d_idv, "cluster_2d_idv[ncluster]/I");
155 tree->Branch("cluster_2d_highest", &cluster_2d_highest,"cluster_2d_highest[ncluster]/I");
156
157}
158
166void TestClusterWithHit::reset_hit() {
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}
177
178void TestClusterWithHit::reset_cluster_1d() {
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}
213
214void TestClusterWithHit::reset_cluster_2d() {
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}
240
241
242
243void TestClusterWithHit::reset_map() {
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}
251
252
253bool TestClusterWithHit::selDigi(CgemDigiCol::iterator iter, int sel)
254{
255 if(sel==0) return true;
256 else {
257 double time = (*iter)->getTime_ns();
258 bool timeIsGood=true;
259 if(time<m_minDigiTime||time>m_maxDigiTime) timeIsGood=false;
260
261 double Q = (*iter)->getCharge_fc();
262 bool chargeIsGood=true;
263 if(Q<myQMin) chargeIsGood=false;
264
265 const Identifier ident = (*iter)->identify();
266 int layer = CgemID::layer(ident);
267 int sheet = CgemID::sheet(ident);
268 int strip = CgemID::strip(ident);
269 int view = 1;
270 bool is_xstrip = CgemID::is_xstrip(ident);
271 if(is_xstrip == true) view = 0;
272 int quality = lutreader->GetQuality(layer, sheet, view, strip);
273 bool qualityIsGood=true;
274 if(quality!=0) qualityIsGood=false;
275
276 if(sel==1) return timeIsGood&&chargeIsGood;
277 else if(sel==-1) return !timeIsGood&&chargeIsGood;
278 else if(sel==2) return timeIsGood&&chargeIsGood&&qualityIsGood;
279 }
280 return false;
281
282}
283
284
285
286
287
288std::map<int, std::vector < int > > * TestClusterWithHit::GetStripToHitMap(int ilayer, int isheet, int iview) {
289
290 if(iview==0) {
291 if(ilayer == 0 && isheet == 0) return &map_L1_S1_stripx_to_hit;
292 else if(ilayer == 1 && isheet == 0) return &map_L2_S2_stripx_to_hit;
293 else if(ilayer == 1 && isheet == 1) return &map_L2_S1_stripx_to_hit;
294 }
295 else if(iview==1) {
296 if(ilayer == 0 && isheet == 0) return &map_L1_S1_stripv_to_hit;
297 else if(ilayer == 1 && isheet == 0) return &map_L2_S2_stripv_to_hit;
298 else if(ilayer == 1 && isheet == 1) return &map_L2_S1_stripv_to_hit;
299 }
300 else return NULL;
301}
302
303StatusCode TestClusterWithHit::execute(){
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}
601
602StatusCode TestClusterWithHit::finalize(){
603 MsgStream log(msgSvc(),name());
604 log << MSG::INFO << "TestClusterWithHit finalize()" << endreq;
605 output->Write();
606 output->Close();
607 return StatusCode::SUCCESS;
608}
609
610
time
Double_t time
Definition DataBase/tau_mode.c:7
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
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
MAXNOFHITS
#define MAXNOFHITS
Definition TestClusterWithHit.h:22
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
CgemLUTReader::GetQuality
int GetQuality(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:378
ICgemGeomSvc::getReadoutPlane
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
TestClusterWithHit::TestClusterWithHit
TestClusterWithHit(const std::string &name, ISvcLocator *pSvcLocator)
Definition TestClusterWithHit.cxx:44