CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
TestTrack.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
32#include "Identifier/CgemID.h"
33
35#include "RawEvent/DigiEvent.h"
38#include "GaudiKernel/SmartDataPtr.h"
39
41
42#include "TMath.h"
43#include "TVector3.h"
44//using namespace std;
45
46TestTrack::TestTrack(const std::string& name, ISvcLocator* pSvcLocator):
47 Algorithm(name,pSvcLocator){
48
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);
55
56
57}
58
60 // delete m_SvcCgem;
61 // delete f;
62 // delete Tdigi;
63 // delete anode;
64 // delete tree;
65 // delete output;
66 delete lutreader;
67}
68
69void TestTrack::DefineHitTree() {
70
71 tree->Branch("event", &event, "event/I");
72 tree->Branch("nhit", &nhit, "nhit/I");
73 // from geometry
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");
79 // from LUT
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");
85 // from signal
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");
90
91
92}
93
94void TestTrack::DefineClusterTree() {
95
96 tree->Branch("ncluster", &ncluster,"ncluster/I");
97 // 1d
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");
124
125 // 2d
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");
142
143}
144
145
146void TestTrack::DefineTrackTree() {
147
148 // parameters
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");
154
155 // clusters
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");
160
161 // layer/sheet under test (if none these are -1)
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");
164
165 // poca
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");
169
170 // intersections
171 // L1
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");
182 // L2
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");
193
194 // Incident Angles
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");
197
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");
200
201
202}
203
204
205
207 reset_map();
210 reset_hit();
211 reset_track();
212}
213
215 // event = 0;
216 nhit = 0;
217
218 for(int ihit = 0; ihit < MAXNOFHITS; ihit++) {
219 hit_strip[ihit] = -1;
220 hit_view[ihit] = -1;
221 hit_layer[ihit] = -1;
222 hit_sheet[ihit] = -1;
223 hit_length[ihit] = -1;
224 hit_channel[ihit] = -1;
225 hit_roc[ihit] = -1;
226 hit_feb[ihit] = -1;
227 hit_tiger[ihit] = -1;
228 hit_chip[ihit] = -1;
229 hit_t[ihit] = 0;
230 hit_q[ihit] = 0;
231 hit_quality[ihit] = 0;
232 hit_selgooddigi[ihit] = 0;
233 }
234}
235
237
238 ncluster_1d = 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;
245
246 for(int iclu = 0; iclu < MAXNOFCLUSTERS; iclu++)
247 {
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;
266 }
267
268
269
270}
271
273
274 ncluster_2d = 0;
275 ncluster_2d_L1_S1 = 0;
276 ncluster_2d_L2_S1 = 0;
277 ncluster_2d_L2_S2 = 0;
278
279 for(int iclu = 0; iclu < MAXNOFCLUSTERS; iclu++)
280 {
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;
296 }
297}
298
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();
306}
307
309
310 // parameters
311 track_dr = -9999;
312 track_phi0 = -9999;
313 track_dz = -9999;
314 track_tanL = -9999;
315 track_chi2 = 9999;
316
317 // clusters
318 track_nfitpoint = 0;
319 for(int iclus=0; iclus < MAXNOFFITPOINT; iclus++) {
320 track_clusterid[iclus] = -1;
321 track_layerid[iclus] = -1;
322 track_sheetid[iclus] = -1;
323 }
324
325 // planes under test
326 track_test_layerid = -1;
327 track_test_sheetid = -1;
328
329 // poca
330 track_xpoca_glo = -9999;
331 track_ypoca_glo = -9999;
332 track_zpoca_glo = -9999;
333
334 // intersections
335 // L1
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;
342 // L2
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;
349
350
351}
352
353bool TestTrack::selDigi(CgemDigiCol::iterator iter, int sel)
354{
355 if(sel==0) return true;
356 else {
357 double time = (*iter)->getTime_ns();
358 bool timeIsGood=true;
359 if(time<m_minDigiTime||time>m_maxDigiTime) timeIsGood=false;
360
361 double Q = (*iter)->getCharge_fc();
362 bool chargeIsGood=true;
363 if(Q<myQMin) chargeIsGood=false;
364
365 const Identifier ident = (*iter)->identify();
366 int layer = CgemID::layer(ident);
367 int sheet = CgemID::sheet(ident);
368 int strip = CgemID::strip(ident);
369 int view = 1;
370 bool is_xstrip = CgemID::is_xstrip(ident);
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;
375
376 if(sel==1) return timeIsGood&&chargeIsGood;
377 else if(sel==-1) return !timeIsGood&&chargeIsGood;
378 else if(sel==2) return timeIsGood&&chargeIsGood&&qualityIsGood;
379 }
380 return false;
381
382}
383
384
385std::map<int, std::vector < int > > * TestTrack::GetStripToHitMap(int ilayer, int isheet, int iview) {
386
387 if(iview==0) {
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;
391 }
392 else if(iview==1) {
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;
396 }
397
398 return NULL;
399}
400
401
403 MsgStream log(msgSvc(), name());
404 log << MSG::INFO << "TestTrack initialize()" << endreq;
405
406
407 // CgemGeomSvc
408 StatusCode sc;
409 sc = service("CgemGeomSvc", m_SvcCgem);
410 if(sc != StatusCode::SUCCESS) {
411 log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
412 return StatusCode::FAILURE;
413 }
414
415 anode = m_SvcCgem->getReadoutPlane(0, 0);
416 anode_mid_gap_L1 = anode->getMidRAtGap();
417 anode = m_SvcCgem->getReadoutPlane(1, 0);
418 anode_mid_gap_L2 = anode->getMidRAtGap();
419
420 midplane = m_SvcCgem->getMidDriftPtr();
421 alignment = m_SvcCgem->getAlignPtr() ;
422
423
424 // LUT
425 lutreader = new CgemLUTReader(lutfile);
426 lutreader->ReadLUT();
427
428 output = new TFile("track.root", "RECREATE");
429 tree = new TTree("tree","hit, cluster and track info tree");
430 DefineHitTree();
431 DefineClusterTree();
432 DefineTrackTree();
433
434 event=0;
435
436 return StatusCode::SUCCESS;
437
438}
439
441
442 MsgStream log(msgSvc(), name());
443
444 // if(event%1000==0) cout << "--------->TestTrack::execute " << event << endl;
445 reset();
446
447 //interface to event data service
448 ISvcLocator* svcLocator = Gaudi::svcLocator();
449 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
450 if (sc.isFailure())
451 cout<<"Could not accesss EventDataSvc!"<<endl;
452 // -------------------------
453 //retrieve CGEM hits from TDS
454 SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
455 if(!aDigiCol)
456 cout<<"Could not retrieve CGEM digi collection"<<endl;
457 nhit = aDigiCol->size();
458 // cout << "total hit " << nhit << endl;
459
460 // loop on hits
461 int ihit=0;
462
463 CgemDigiCol::iterator iDigiCol;
464 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
465 {
466 const Identifier ident = (*iDigiCol)->identify();
467 int strip = CgemID::strip(ident);
468 int layer = CgemID::layer(ident);
469 int sheet = CgemID::sheet(ident);
470 int view = 1;
471 bool is_xstrip = CgemID::is_xstrip(ident);
472 if(is_xstrip == true) view = 0;
473 double charge = (*iDigiCol)->getCharge_fc();
474 double time = (*iDigiCol)->getTime_ns();
475
476 hit_strip[ihit] = strip;
477 hit_layer[ihit] = layer;
478 hit_sheet[ihit] = sheet;
479 hit_view[ihit] = view;
480 hit_t[ihit] = time;
481 hit_q[ihit] = charge;
482
483 anode = m_SvcCgem->getReadoutPlane(layer, sheet);
484 if(view == 1) hit_length[ihit] = anode->getVStripLength(strip);
485 else hit_length[ihit] = anode->getLength();
486
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);
493
494 hit_channel[ihit] = channel;
495 hit_roc[ihit] = roc;
496 hit_feb[ihit] = feb;
497 hit_tiger[ihit] = tiger;
498 hit_chip[ihit] = chip;
499 hit_quality[ihit] = quality;
500
501 if(selDigi(iDigiCol, m_selGoodDigi)==true) {
502 hit_selgooddigi[ihit] = 1;
503
504 // cout << "HIT " << ihit << " strip " << strip << " layer " << layer << " sheet " << sheet << " view " << view << endl;
505 // map strips to hits (maybe more hits for the same strip)
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);
510 }
511 else {
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);
516 }
517 }
518 // --------------------------------------------------------
519 ihit++;
520 }
521
522
523
524 // -------------------------------------------------------------------
525 //retrieve CGEM clusters from TDS
526 SmartDataPtr<RecCgemClusterCol> aClusterCol(m_evtSvc,"/Event/Recon/RecCgemClusterCol");
527 if(!aClusterCol)
528 cout<<"Could not retrieve CGEM cluster collection"<<endl;
529 int nclu = aClusterCol->size();
530 //cout << "total cluster " << nclu << endl;
531 if(nclu > MAXNOFCLUSTERS) {
532 std::cout << "nclu " << nclu << std::endl;
533 nclu = MAXNOFCLUSTERS; // CHECK HACK
534 }
535
536 int nclusterx = 0;
537 int nclusterv = 0;
538 int nclusterxv = 0;
539
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.;
546
547 ncluster = 0;
548
549 RecCgemClusterCol::iterator iClusterCol;
550 for(iClusterCol=aClusterCol->begin(); iClusterCol!=aClusterCol->end(); iClusterCol++)
551 {
552
553 if(ncluster==MAXNOFCLUSTERS) break;
554
555 int flag = (*iClusterCol)->getflag(); // 0=x 1=v 2=xv
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();
570 double a_tpc = 0; //(*iClusterCol)->getparamA_mTPC(); CHECK
571 double b_tpc = 0; //(*iClusterCol)->getparamB_mTPC(); CHECK
572
573 if(sheetid == -1) {
574 ncluster++;
575 continue; // CHECK HACK FOR NOW
576 }
577
578 // anode mid gap radius
579 anode = m_SvcCgem->getReadoutPlane(layerid, sheetid);
580 double anode_mid_gap = anode->getMidRAtGap();
581
582 // beginning and ending flags
583 // if 2D cluster --> they are the x & v 1D clusters
584 // if 1D cluster --> they are the first and last contigous strips
585 int flagB = (*iClusterCol)->getclusterflagb();
586 int flagE = (*iClusterCol)->getclusterflage();
587 // ------------------------------------------------------------------- 1D cluster
588 // positions CC & uTPC
589 if(flag==0 || flag == 1) {
590 cluster_1d_ID[ncluster] = clusterid;
591 // time & charge
592 cluster_1d_t[ncluster] = 0; // CHECK
593 cluster_1d_q[ncluster] = edep;
594 // layer/sheet/view
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;
607 // from strip1 to strip2
608 cluster_1d_strip1[ncluster] = flagB;
609 cluster_1d_strip2[ncluster] = flagE;
610 int size = flagE - flagB + 1;
611 cluster_1d_size[ncluster] = size;
612
613 // get the hits
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;
616 /**
617 cout << "1D_CLUSTER index " << clusterid << " or " << ncluster
618 << " layer " << layerid << " sheet " << sheetid
619 << " view " << cluster_1d_view[ncluster] << endl;
620 cout << " from " << flagB << " to " << flagE << endl;
621 cout << "hits: ";
622 **/
623
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]; // take only the first hit CHECK THIS HAS TO BE CHANGED!
630 }
631 // -----------------
632 ncluster_1d++;
633 }
634 // ----------
635 else if(flag==2 || flag==3) { // ------------------------------------- 2D - xv view
636 cluster_2d_ID[ncluster] = clusterid;
637 // time & charge
638 cluster_2d_t[ncluster] = 0; // CHECK
639 cluster_2d_q[ncluster] = edep;
640 // layer/sheet/view
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;
652 // composing 2D_clusters
653 cluster_2d_idx[ncluster] = flagB;
654 cluster_2d_idv[ncluster] = flagE;
655
656 /** cout << "2D_CLUSTER index " << clusterid << " or " << ncluster
657 << " layer " << layerid << " sheet " << sheetid
658 << " view " << cluster_2d_view[ncluster] << endl;
659 cout << " idx " << flagB << " idv " << flagE << endl;
660 **/
661
662 if(flag==2) { // HACK THIS
663 // find the highest charge cluster
664 if(layerid==0) {
665 if(edep > tmp_charge_L1_S1) {
666 tmp_charge_L1_S1 = edep; tmp_cluster_L1_S1 = ncluster;
667 }
668 }
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;
672 }
673 }
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;
677 }
678 }
679 }
680 // -----------
681 ncluster_2d++;
682 }
683 // cluster counters
684 ncluster++;
685 if(flag == 0) {
686 nclusterx++;
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++;
691 }
692 }
693 else if(flag == 1) {
694 nclusterv++;
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++;
699 }
700 }
701 else if(flag == 2) {
702 nclusterxv++;
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++;
707 }
708 }
709 // --------
710 }
711
712 // std::cout << "highest charge cluster " << tmp_cluster_L1_S1 << " " << tmp_cluster_L2_S1 << " " << tmp_cluster_L2_S2 << std::endl;
713 // std::cout << "charge " << tmp_charge_L1_S1 << " " << tmp_charge_L2_S1 << " " << tmp_charge_L2_S2 << std::endl;
714
715 // set the highest charge 2D - cluster
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;
719
720 // -------------------------------------------------------------------
721 //retrieve CGEM tracks from TDS
722 SmartDataPtr<RecMdcTrackCol> aTrackCol(m_evtSvc,"/Event/Recon/RecMdcTrackCol");
723 if(!aTrackCol) {
724 // cout<<"Could not retrieve CGEM track collection"<<endl;
725 }
726 else {
727 // cout << "total track " << aTrackCol->size() << endl;
728
729 RecMdcTrackCol::iterator iTrackCol;
730 for(iTrackCol=aTrackCol->begin(); iTrackCol!=aTrackCol->end(); iTrackCol++)
731 {
732 RecMdcTrack *track = *iTrackCol;
733
734 /**
735 cout << "track id " << track->trackId() << endl;
736 cout << "dr " << track->helix(0) << endl;
737 cout << "phi0 " << track->helix(1) << endl;
738 cout << "dz " << track->helix(3) << endl;
739 cout << "tanL " << track->helix(4) << endl;
740 **/
741
742 if(track->trackId() != 0) continue;
743
744 track_chi2 = track->chi2();
745 // cout << "chi2 " << track_chi2 << endl;
746 track_dr = track->helix(0);
747 track_phi0 = track->helix(1);
748 track_dz = track->helix(3);
749 track_tanL = track->helix(4);
750
751 ClusterRefVec vecclusters = track->getVecClusters();
752 track_nfitpoint = vecclusters.size(); // nof 2D clusters, NOTE! track->getNcluster() gives you the number of 1D clusters
753
754 // is the plane under test?
755 bool istestplane[2][2];
756 for(int ilay=0; ilay<MAXNOFLAYER; ilay++) for(int ishe=0; ishe<MAXNOFSHEET; ishe++) istestplane[ilay][ishe] = true;
757
758 for(int iclus=0; iclus < track_nfitpoint; iclus++) {
759 RecCgemCluster *tcluster = vecclusters.at(iclus);
760 double phi = tcluster->getrecphi();
761 double z = tcluster->getRecZ();
762 track_clusterid[iclus] = tcluster->getclusterid();
763
764 track_layerid[iclus] = tcluster->getlayerid();
765 track_sheetid[iclus] = tcluster->getsheetid();
766 if(phi > 0) track_sheetid[iclus] = 1;
767
768 // cout << "iclus " << iclus << " track_clusterid " << track_clusterid[iclus] << " phi " << phi
769 // << " layer " << track_layerid[iclus] << " sheet " << track_sheetid[iclus] << endl;
770
771 istestplane[track_layerid[iclus]][track_sheetid[iclus]] = false;
772 }
773
774 for(int ilay=0; ilay<MAXNOFLAYER; ilay++) {
775 for(int ishe=0; ishe<MAXNOFSHEET; ishe++) {
776
777 // cout << "situazione ilay " << ilay << " ishe " << ishe << " istestplane " << istestplane[ilay][ishe] << endl;
778
779 if(istestplane[ilay][ishe] == true) {
780 track_test_layerid = ilay;
781 track_test_sheetid = ishe;
782 }
783 }
784 }
785
786 // cout << "UNDER TEST layer " << track_test_layerid << " sheet " << track_test_sheetid << endl;
787
788 // intersections - poca and intersections are computed in the aligned frame
789 double xP; double yP; double zP;
790 ComputePOCA(xP, yP, zP);
791
792 // L1
793 double xP1; double yP1; double zP1; double phiP1; double vP1;
794 double xP2; double yP2; double zP2; double phiP2; double vP2;
795
796 double angCR_xy_L1, angCR_yz_L1;
797
798 bool gotit1 = ComputeIntersection(0, xP1, yP1, zP1, phiP1, vP1, xP2, yP2, zP2, phiP2, vP2, angCR_xy_L1, angCR_yz_L1); //1->up 2->down
799
800 ang_xy_L1 = angCR_xy_L1;
801 ang_yz_L1 = angCR_yz_L1;
802
803 // L2
804 double xP3; double yP3; double zP3; double phiP3; double vP3;
805 double xP4; double yP4; double zP4; double phiP4; double vP4;
806
807 double angCR_xy_L2, angCR_yz_L2;
808
809 bool gotit2 = ComputeIntersection(1, xP3, yP3, zP3, phiP3, vP3, xP4, yP4, zP4, phiP4, vP4, angCR_xy_L2, angCR_yz_L2); //1->up 2->down
810
811 ang_xy_L2 = angCR_xy_L2;
812 ang_yz_L2 = angCR_yz_L2;
813
814 if(gotit1==false || gotit2==false) cout << "TestTrack: error in intersection calculation. Intersection on L1 " << gotit1 << ", on L2 " << gotit2 << endl;
815
816 // output
817 // poca
818 track_xpoca_glo = xP; track_ypoca_glo = yP; track_zpoca_glo = zP;
819 // intersections
820 // L1
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;
823 // L2
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;
826 }
827
828
829 }
830 tree->Fill();
831
832 event++;
833
834 return StatusCode::SUCCESS;
835}
836
838 MsgStream log(msgSvc(),name());
839 log << MSG::INFO << "TestTrack finalize()" << endreq;
840 output->Write();
841 output->Close();
842 return StatusCode::SUCCESS;
843}
844
845// computation of P.O.C.A. = Point Of Closest Approach of the found
846// track to the origin (0, 0, 0,) in the global, aligned reference frame
847// --> xP, yP, zP is in the aligned frame!
848void TestTrack::ComputePOCA(double &xP, double &yP, double &zP)
849{
850 // poca
851 xP = track_dr * cos(track_phi0);
852 yP = track_dr * sin(track_phi0);
853 zP = track_dz;
854}
855
856// Compute the intersection of the fitted track on the layerid-ith layer
857// * layerid = 0, 1, 2 means L1, L2, L3
858// * if align flag is on --> the intersection x, y, z is in the alighed frame
859// otherwise the intersection x, y, z are in the non aligned frame
860// * phi1, v1 are always in the local frame of the layer
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)
865{
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];
869 bool gotit;
870
871 if(align_flag==true) {
872 // sheet bottom
873 gotit = midplane->getPointAligned(layerid, 0, linefit, posup_temp, posdown, phivup_temp, phivdown);
874 // sheet top
875 gotit = midplane->getPointAligned(layerid, 1, linefit, posup, posdown_temp, phivup, phivdown_temp);
876 }
877 else gotit = midplane->getPointIdealGeom(layerid, linefit, posup, posdown, phivup, phivdown);
878
879 x1 = posup.x();
880 y1 = posup.y();
881 z1 = posup.z();
882 phi1 = phivup[0];
883 v1 = phivup[1];
884
885 x2 = posdown.x();
886 y2 = posdown.y();
887 z2 = posdown.z();
888 phi2 = phivdown[0];
889 v2 = phivdown[1];
890
891 Hep3Vector CosmicRay_xy(x1-x2, y1-y2, 0);
892 Hep3Vector CosmicRay_yz( 0, y1-y2, z1-z2);
893
894 Hep3Vector NormalToCRayInPoint(x1, y1, 0);
895 Hep3Vector zAxis(0, 0, 1);
896
897 angCR_xy = CosmicRay_xy.angle(NormalToCRayInPoint);
898 angCR_yz = CosmicRay_yz.angle(zAxis);
899
900 /**
901 cout << "layerid " << layerid << endl;
902 cout << "up phi "<< phivup[0] << " da y/x " << TMath::ATan2(posup.y(), posup.x()) << endl;
903 cout << "down phi "<< phivdown[0] << " da y/x " << TMath::ATan2(posdown.y(), posdown.x()) << endl;
904 **/
905
906 return gotit;
907
908}
909
910
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
Double_t phi2
Double_t phi1
Double_t time
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
**********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
SmartRefVector< RecCgemCluster > ClusterRefVec
Definition: RecMdcTrack.h:28
IMessageSvc * msgSvc()
#define MAXNOFSHEET
Definition: TestTrack.h:26
#define MAXNOFFITPOINT
Definition: TestTrack.h:24
#define MAXNOFCLUSTERS
Definition: TestTrack.h:22
#define MAXNOFHITS
Definition: TestTrack.h:23
#define MAXNOFLAYER
Definition: TestTrack.h:25
double getVStripLength(int V_ID) const
static int strip(const Identifier &id)
Definition: CgemID.cxx:83
static int sheet(const Identifier &id)
Definition: CgemID.cxx:77
static int layer(const Identifier &id)
Definition: CgemID.cxx:71
static bool is_xstrip(const Identifier &id)
Definition: CgemID.cxx:64
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
Definition: DstMdcTrack.h:66
const int trackId() const
Definition: DstMdcTrack.h:52
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
Definition: RecMdcTrack.h:70
void reset_cluster_2d()
Definition: TestTrack.cxx:272
void reset_cluster_1d()
Definition: TestTrack.cxx:236
StatusCode finalize()
Definition: TestTrack.cxx:837
void reset()
Definition: TestTrack.cxx:206
TestTrack(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TestTrack.cxx:46
void reset_map()
Definition: TestTrack.cxx:299
void reset_track()
Definition: TestTrack.cxx:308
StatusCode initialize()
Definition: TestTrack.cxx:402
StatusCode execute()
Definition: TestTrack.cxx:440
void reset_hit()
Definition: TestTrack.cxx:214