CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
TestHit.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 "Identifier/CgemID.h"
31
32#include "RawEvent/RawDataUtil.h"
33#include "RawEvent/DigiEvent.h"
34#include "ReconEvent/ReconEvent.h"
35#include "EventModel/EventHeader.h"
36#include "GaudiKernel/SmartDataPtr.h"
37
38#include "ReadCosmicRayData/TestHit.h"
39#include "TMath.h"
40//using namespace std;
41
42TestHit::TestHit(const std::string& name, ISvcLocator* pSvcLocator):
43 Algorithm(name,pSvcLocator){
44
45 declareProperty("LUTfile", lutfile = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_from_17_to_17.root");
46 declareProperty("TimeMin", timemin = -8900); // noise or signal
47 declareProperty("TimeMax", timemax = -8500); //
48
49}
50
51TestHit::~TestHit(){
52
53 // delete tree;
54 // delete output;
55 // delete lutreader;
56}
57
58StatusCode TestHit::initialize(){
59
60 MsgStream log(msgSvc(), name());
61 log << MSG::INFO << "TestHit initialize()" << endreq;
62
63 // CgemGeomSvc
64 StatusCode sc;
65
66 sc = service("CgemGeomSvc", m_geoSvc);
67 if(sc != StatusCode::SUCCESS) {
68 log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
69 return StatusCode::FAILURE;
70 }
71
72 event = 0;
73
74 output = new TFile("hit.root", "RECREATE");
75 tree = new TTree("tree","hit info tree");
76 tree->Branch("event", &event, "event/I");
77 tree->Branch("nhit", &nhit, "nhit/I");
78 tree->Branch("nhit_L1_S1_x", &nhit_L1_S1_x,"nhit_L1_S1_x/I");
79 tree->Branch("nhit_L2_S1_x", &nhit_L2_S1_x,"nhit_L2_S1_x/I");
80 tree->Branch("nhit_L2_S2_x", &nhit_L2_S2_x,"nhit_L2_S2_x/I");
81 tree->Branch("nhit_L1_S1_v", &nhit_L1_S1_v,"nhit_L1_S1_v/I");
82 tree->Branch("nhit_L2_S1_v", &nhit_L2_S1_v,"nhit_L2_S1_v/I");
83 tree->Branch("nhit_L2_S2_v", &nhit_L2_S2_v,"nhit_L2_S2_v/I");
84 //
85 tree->Branch("ntwin_L1_S1_x", &ntwin_L1_S1_x,"ntwin_L1_S1_x/I");
86 tree->Branch("ntwin_L2_S1_x", &ntwin_L2_S1_x,"ntwin_L2_S1_x/I");
87 tree->Branch("ntwin_L2_S2_x", &ntwin_L2_S2_x,"ntwin_L2_S2_x/I");
88 tree->Branch("ntwin_L1_S1_v", &ntwin_L1_S1_v,"ntwin_L1_S1_v/I");
89 tree->Branch("ntwin_L2_S1_v", &ntwin_L2_S1_v,"ntwin_L2_S1_v/I");
90 tree->Branch("ntwin_L2_S2_v", &ntwin_L2_S2_v,"ntwin_L2_S2_v/I");
91 // from geometry
92 tree->Branch("hit_strip", &hit_strip, "hit_strip[nhit]/I");
93 tree->Branch("hit_view", &hit_view, "hit_view[nhit]/I");
94 tree->Branch("hit_layer", &hit_layer, "hit_layer[nhit]/I");
95 tree->Branch("hit_sheet", &hit_sheet, "hit_sheet[nhit]/I");
96 tree->Branch("hit_length", &hit_length, "hit_length[nhit]/D");
97 // from LUT
98 tree->Branch("hit_channel", &hit_channel,"hit_channel[nhit]/I");
99 tree->Branch("hit_roc", &hit_roc, "hit_roc[nhit]/I");
100 tree->Branch("hit_feb", &hit_feb, "hit_feb[nhit]/I");
101 tree->Branch("hit_tiger", &hit_tiger, "hit_tiger[nhit]/I");
102 tree->Branch("hit_chip", &hit_chip, "hit_chip[nhit]/I");
103 // from signal
104 tree->Branch("hit_t", &hit_t, "hit_t[nhit]/D");
105 tree->Branch("hit_q", &hit_q, "hit_q[nhit]/D");
106 tree->Branch("hit_saturated", &hit_saturated, "hit_saturated[nhit]/I");
107 tree->Branch("hit_quality", &hit_quality, "hit_quality[nhit]/I");
108
109 // LUT
110 lutreader = new CgemLUTReader(lutfile);
111 lutreader->ReadLUT();
112
113 return StatusCode::SUCCESS;
114}
115
116void TestHit::reset() {
117
118 event = 0;
119 nhit = 0;
120 nhit_L1_S1_x= 0 ;
121 nhit_L2_S1_x= 0 ;
122 nhit_L2_S2_x= 0 ;
123 nhit_L1_S1_v= 0 ;
124 nhit_L2_S1_v= 0 ;
125 nhit_L2_S2_v= 0 ;
126
127 ntwin_L1_S1_x= 0 ;
128 ntwin_L2_S1_x= 0 ;
129 ntwin_L2_S2_x= 0 ;
130 ntwin_L1_S1_v= 0 ;
131 ntwin_L2_S1_v= 0 ;
132 ntwin_L2_S2_v= 0 ;
133
134 for(int ihit = 0; ihit < MAXNOFHIT; ihit++) {
135 hit_strip[ihit] = -1;
136 hit_view[ihit] = -1;
137 hit_layer[ihit] = -1;
138 hit_sheet[ihit] = -1;
139 hit_length[ihit] = 0.;
140
141 hit_channel[ihit] = -1;
142 hit_roc[ihit] = -1;
143 hit_feb[ihit] = -1;
144 hit_tiger[ihit] = -1;
145 hit_chip[ihit] = -1;
146 hit_quality[ihit] = -1;
147
148 hit_t[ihit] = 0 ;
149 hit_q[ihit] = 0 ;
150 hit_saturated[ihit] = 0;
151 }
152
153}
154
155StatusCode TestHit::execute(){
156
157 MsgStream log(msgSvc(), name());
158 // cout << "->TestHit::execute" << endl;
159
160 //interface to event data service
161 ISvcLocator* svcLocator = Gaudi::svcLocator();
162 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
163 if (sc.isFailure())
164 cout<<"Could not accesss EventDataSvc!"<<endl;
165
166 // -------------------
167
168
169 // CgemGeoReadoutPlane* anode = m_SvcCgem->getReadoutPlane(0, 0);
170 // anode_radius_L1_x = anode->getRX();
171 // anode_radius_L1_v = anode->getRV();
172 // anode_mid_gap_L1 = anode->getMidRAtGap();
173 // anode = m_SvcCgem->getReadoutPlane(1, 0);
174 // anode_radius_L2_x = anode->getRX();
175 // anode_radius_L2_v = anode->getRV();
176 // anode_mid_gap_L2 = anode->getMidRAtGap();
177
178 // -------------------
179
180
181 reset();
182 // -------------
183
184 // retrieve CGEM digits from TDS
185 SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
186 if(!aDigiCol)
187 cout<<"Could not retrieve CGEM digi collection"<<endl;
188
189 nhit = aDigiCol->size();
190 CgemDigiCol::iterator iDigiCol;
191
192 // event = m_evtSvc->GetEventID(); // CHECK
193
194 // loop on hits
195 int ihit=0;
196 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
197 {
198 const Identifier ident = (*iDigiCol)->identify();
199
200 int layer = CgemID::layer(ident);
201 int sheet = CgemID::sheet(ident);
202 int strip = CgemID::strip(ident);
203 int view = 1;
204 bool is_xstrip = CgemID::is_xstrip(ident);
205 if(is_xstrip == true) view = 0;
206 double charge = (*iDigiCol)->getCharge_fc();
207 double time = (*iDigiCol)->getTime_ns();
208
209 hit_strip[ihit] = strip;
210 hit_layer[ihit] = layer;
211 hit_sheet[ihit] = sheet;
212 hit_view[ihit] = view;
213
214 anode = m_geoSvc->getReadoutPlane(layer, sheet);
215 if(view == 1) hit_length[ihit] = anode->getVStripLength(strip);
216 else hit_length[ihit] = anode->getLength();
217
218 hit_channel[ihit] = lutreader->GetChannel(layer, sheet, view, strip);
219 hit_roc[ihit] = lutreader->GetROC(layer, sheet, view, strip);
220 hit_feb[ihit] = lutreader->GetFEB(layer, sheet, view, strip);
221 hit_tiger[ihit] = lutreader->GetTIGER(layer, sheet, view, strip);
222 hit_chip[ihit] = lutreader->GetChip(layer, sheet, view, strip);
223 hit_quality[ihit] = lutreader->GetQuality(layer, sheet, view, strip);
224
225 // std::cout << ihit << " strip " << strip << " ch " << hit_channel[ihit] << " roc " << hit_roc[ihit] << " feb " << hit_feb[ihit] << " tg " << hit_tiger[ihit] << std::endl;
226
227 hit_t[ihit] = time;
228 hit_q[ihit] = charge;
229 hit_saturated[ihit] = 0; // CHECK
230
231 if(is_xstrip == true) {
232 if(layer==0 && sheet==0) nhit_L1_S1_x++;
233 else if(layer==1 && sheet==0) nhit_L2_S1_x++;
234 else if(layer==1 && sheet==1) nhit_L2_S2_x++;
235 }
236 else {
237 if(layer==0 && sheet==0) nhit_L1_S1_v++;
238 else if(layer==1 && sheet==0) nhit_L2_S1_v++;
239 else if(layer==1 && sheet==1) nhit_L2_S2_v++;
240 }
241
242 // --- time window
243 if(time > timemin && time < timemax) {
244 if(is_xstrip == true) {
245 if(layer==0 && sheet==0) ntwin_L1_S1_x++;
246 else if(layer==1 && sheet==0) ntwin_L2_S1_x++;
247 else if(layer==1 && sheet==1) ntwin_L2_S2_x++;
248 }
249 else {
250 if(layer==0 && sheet==0) ntwin_L1_S1_v++;
251 else if(layer==1 && sheet==0) ntwin_L2_S1_v++;
252 else if(layer==1 && sheet==1) ntwin_L2_S2_v++;
253 }
254 }
255 // -----------
256 ihit++;
257 }
258
259 tree->Fill();
260 event++;
261
262 return StatusCode::SUCCESS;
263}
264
265StatusCode TestHit::finalize(){
266
267 MsgStream log(msgSvc(),name());
268 log << MSG::INFO << "TestHit finalize()" << endreq;
269 output->Write();
270 output->Close();
271
272 return StatusCode::SUCCESS;
273}
274
275
time
Double_t time
Definition: DataBase/tau_mode.c:7
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
MAXNOFHIT
#define MAXNOFHIT
Definition: TestHit.h:19
CgemGeoReadoutPlane::getVStripLength
double getVStripLength(int V_ID) const
Definition: CgemGeoReadoutPlane.h:83
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::GetChip
int GetChip(int ilayer, int isheet, int iview, int istrip)
Definition: CgemLUTReader.cxx:191
CgemLUTReader::GetQuality
int GetQuality(int ilayer, int isheet, int iview, int istrip)
Definition: CgemLUTReader.cxx:371
CgemLUTReader::GetROC
int GetROC(int ilayer, int isheet, int iview, int istrip)
Definition: CgemLUTReader.cxx:173
CgemLUTReader::GetChannel
int GetChannel(int ilayer, int isheet, int iview, int istrip)
Definition: CgemLUTReader.cxx:197
CgemLUTReader::GetTIGER
int GetTIGER(int ilayer, int isheet, int iview, int istrip)
Definition: CgemLUTReader.cxx:179
CgemLUTReader::GetFEB
int GetFEB(int ilayer, int isheet, int iview, int istrip)
Definition: CgemLUTReader.cxx:185
ICgemGeomSvc::getReadoutPlane
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
TestHit::initialize
StatusCode initialize()
Definition: TestHit.cxx:58
TestHit::~TestHit
~TestHit()
Definition: TestHit.cxx:51
TestHit::execute
StatusCode execute()
Definition: TestHit.cxx:155
TestHit::TestHit
TestHit(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TestHit.cxx:42
TestHit::finalize
StatusCode finalize()
Definition: TestHit.cxx:265
TestHit::reset
void reset()
Definition: TestHit.cxx:116