CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
bak_ReadCosmicRayData-00-00-04/src/ReadCosmicRayData.cxx
Go to the documentation of this file.
1/**************************************************************************
2 * CgemBOSS (BESIII Offline Software System) *
3 * *
4 * Author: The BESII Collaboration *
5 * Contributors: Aiqiang Guo, Liangliang Wang *
6 * Date: Dec 10 2018 *
7 * *
8 **************************************************************************/
9/**************************************************************************
10 * Note for 01
11 * Change the reading method for vectorized data
12 * By: Aiqiang Guo
13 * Date: Dec 17 2018
14 **************************************************************************/
15/**************************************************************************
16 * Note for 02
17 * Add CC and mTPC for cluster object
18 * By: Aiqiang Guo
19 * Date: Dec 27 2018
20 **************************************************************************/
21
22// include system lib
23#include <iostream>
24#include <iomanip>
25#include <string>
26#include <cmath>
27
28// Include files
29#include "GaudiKernel/AlgFactory.h"
30#include "GaudiKernel/DataObject.h"
31#include "GaudiKernel/IEventProcessor.h"
32
33#include "GaudiKernel/Incident.h"
34#include "GaudiKernel/IIncidentSvc.h"
35#include "GaudiKernel/Memory.h"
36
37#include <csignal>
38
39// for save digit & cluster
40#include "GaudiKernel/ISvcLocator.h"
41#include "GaudiKernel/IDataProviderSvc.h"
42#include "GaudiKernel/Bootstrap.h"
43#include "GaudiKernel/RegistryEntry.h"
44#include "GaudiKernel/MsgStream.h"
45
46#include "CgemRawEvent/CgemDigi.h"
47#include "CgemRecEvent/RecCgemCluster.h"
48
49#include "Identifier/CgemID.h"
50
51#include "RawEvent/RawDataUtil.h"
52#include "RawEvent/DigiEvent.h"
53#include "ReconEvent/ReconEvent.h"
54#include "EventModel/EventHeader.h"
55#include "GaudiKernel/SmartDataPtr.h"
56
57#include "ReadCosmicRayData/ReadCosmicRayData.h"
58
59//using namespace std;
60ReadCosmicRayData::ReadCosmicRayData(const std::string& name, ISvcLocator* pSvcLocator):
61 Algorithm(name,pSvcLocator){
62
63 declareProperty("Dir_file", Dir_file = "Cosmic_data_01.root");
64 declareProperty("TreeDigi", TreeDigi = "t1");
65 declareProperty("TreeCluster", TreeCluster = "t1");
66 declareProperty("ReadDigi", ReadDigi = true);
67 declareProperty("ReadCluster", ReadCluster = true);
68 declareProperty("DigiSheetID", DigiSheetID = 0);
69 declareProperty("Cut_on_tpc", Cut_on_tpc = false);
70 declareProperty("ClusterSheetID", ClusterSheetID = 0);
71 declareProperty("ClusterRecZ", ClusterRecZ = 0);
72 declareProperty("R_Cluster", R_Cluster = 1.0);
73 declareProperty("Shift_DigitLayerID", Shift_DigitLayerID = 0);
74 declareProperty("Shift_DigitSheetID", Shift_DigitSheetID = 0);
75 declareProperty("Shift_DigitXStripID", Shift_DigitXStripID = 0);
76 declareProperty("Shift_DigitVStripID", Shift_DigitVStripID = 0);
77 declareProperty("Shift_ClusterLayerID", Shift_ClusterLayerID = 0);
78 declareProperty("Shift_ClusterSheetID", Shift_ClusterSheetID = 0);
79 declareProperty("Shift_RecPhi", Shift_RecPhi = 0);
80 declareProperty("Shift_RecV", Shift_RecV = 0);
81 declareProperty("Shift_RecZ", Shift_RecZ = 0);
82 declareProperty("CosmicRayDataSetID", CosmicRayDataSetID = "CR201909");
83
84}
85
87}
88
90 MsgStream log(msgSvc(), name());
91 log << MSG::INFO << "ReadCosmicRayData initialize()" << endreq;
92
93 TString TDir_file(Dir_file);
94 f = new TFile(TDir_file);
95 if(ReadDigi)
96 {
97 TString TTreeDigi(TreeDigi);
98 Tdigi = (TTree*)f->Get(TTreeDigi);
99
100 // Get Cgem digi tree
101 Tdigi->SetBranchAddress("Event", &m_Event_D); // event ID
102 Tdigi->SetBranchAddress("nGemHit", &m_nGemHit); // nof GEM hits
103 // Tdigi->SetBranchAddress("GemHit_nHit", &m_nGemHit); // nof GEM hits it is the same thing as before
104
105 // information on the IDs
106 Tdigi->SetBranchAddress("GemHit_channel", m_channel); // channel no. [0, 63]
107 Tdigi->SetBranchAddress("GemHit_ROC", m_ROC); // ROC no.
108 Tdigi->SetBranchAddress("GemHit_chip", m_chip); // chip no.
109 Tdigi->SetBranchAddress("GemHit_FEB", m_FEB); // FEB no.
110 Tdigi->SetBranchAddress("GemHit_plane", m_plane); // plane
111 Tdigi->SetBranchAddress("GemHit_view", m_view); // view (axial or stereo strips)
112 Tdigi->SetBranchAddress("GemHit_strip", m_strip); // strip no.
113
114 // physical information
115 Tdigi->SetBranchAddress("GemHit_saturated", m_saturated); // is the ASIC channel saturated
116 Tdigi->SetBranchAddress("GemHit_q", m_charge); // charge (fC)
117 Tdigi->SetBranchAddress("GemHit_time", m_time); // time (ns)
118 // Tdigi->SetBranchAddress("GemHit_is_tpc", m_GemHit_is_tpc);
119
120 No_Entries_D = Tdigi->GetEntries();
121 Ind_Entry_D = 0;
122 }
123
124 /**
125 if(ReadCluster)
126 {
127 //Get Cgem cluster tree
128 TString TTreeCluster(TreeCluster);
129 Tcluster = (TTree*)f->Get(TTreeCluster);
130
131 Tcluster->SetBranchAddress("Event", &m_Event_C);
132 Tcluster->SetBranchAddress("GemCluster1d_nCluster", &m_nGemCluster);
133 Tcluster->SetBranchAddress("GemCluster1d_nHit", m_ClusternHit);
134 Tcluster->SetBranchAddress("GemCluster1d_HitIndex", m_ClusterHitIndex);
135 Tcluster->SetBranchAddress("GemCluster1d_plane", m_ClusterLayerID);
136 Tcluster->SetBranchAddress("GemCluster1d_view", m_Flag);
137 Tcluster->SetBranchAddress("GemCluster1d_q", m_EnergyDeposit);
138 Tcluster->SetBranchAddress("GemCluster1d_x", m_Cluster_x);
139 Tcluster->SetBranchAddress("GemCluster1d_z", m_Cluster_z);
140 Tcluster->SetBranchAddress("GemCluster1d_x_cc", m_Cluster_x_cc);
141 Tcluster->SetBranchAddress("GemCluster1d_x_tpc", m_Cluster_x_tpc);
142 Tcluster->SetBranchAddress("GemCluster1d_z_cc", m_Cluster_z_cc);
143 Tcluster->SetBranchAddress("GemCluster1d_z_tpc", m_Cluster_z_tpc);
144 //Tcluster->SetBranchAddress("ClusterSheetID", m_ClusterSheetID);
145 //Tcluster->SetBranchAddress("ClusterFlagB", m_ClusterFlagB);
146 //Tcluster->SetBranchAddress("ClusterFlagE", m_ClusterFlagE);
147 //Tcluster->SetBranchAddress("RecV", m_RecV);
148 //Tcluster->SetBranchAddress("RecZ", m_RecZ);
149
150 No_Entries_C = Tcluster->GetEntries();
151 Ind_Entry_C = 0;
152 }
153 **/
154
155 return StatusCode::SUCCESS;
156}
157
158
159int ReadCosmicRayData::TranslateDigitLayerID(int Input_LayerID)
160{
161 int ShiftValue = Shift_DigitLayerID;
162 return Input_LayerID+ShiftValue;
163}
164
165int ReadCosmicRayData::TranslateDigitSheetID(int Input_SheetID)
166{
167 int ShiftValue = Shift_DigitSheetID;
168 return Input_SheetID+ShiftValue;
169}
170
171int ReadCosmicRayData::TranslateDigitXStripID(int Input_StripID)
172{
173 int ShiftValue = Shift_DigitXStripID;
174 return Input_StripID+ShiftValue;
175}
176
177int ReadCosmicRayData::TranslateDigitVStripID(int Input_StripID)
178{
179 int ShiftValue = Shift_DigitVStripID;
180 return Input_StripID+ShiftValue;
181}
182
183int ReadCosmicRayData::TranslateDigitStripID(int Input_StripID, int StripType)
184{
185 int Output_StripID = -1;
186 if(StripType==0) Output_StripID = TranslateDigitXStripID(Input_StripID);
187 if(StripType==1) Output_StripID = TranslateDigitVStripID(Input_StripID);
188 return Output_StripID;
189}
190
191int ReadCosmicRayData::TranslateDigitStripType(int Input_StripType)
192{
193 return Input_StripType;
194}
195
196int ReadCosmicRayData::TranslateClusterLayerID(int Input_LayerID)
197{
198 int ShiftValue = Shift_ClusterLayerID;
199 return Input_LayerID+ShiftValue;
200}
201
202int ReadCosmicRayData::TranslateClusterSheetID(int Input_SheetID)
203{
204 int ShiftValue = Shift_ClusterSheetID;
205 return Input_SheetID+ShiftValue;
206}
207
208int ReadCosmicRayData::TranslateClusterFlag(int Input_Flag)
209{
210 return Input_Flag-2;
211}
212
213double ReadCosmicRayData::TranslateRecPhi(double Input_RecPhi)
214{
215 double ShiftValue = Shift_RecPhi;
216 return Input_RecPhi+ShiftValue;
217}
218
219double ReadCosmicRayData::TranslateRecV(double Input_RecV)
220{
221 double ShiftValue = Shift_RecV;
222 return Input_RecV+ShiftValue;
223}
224
225double ReadCosmicRayData::TranslateRecZ(double Input_RecZ)
226{
227 double ShiftValue = Shift_RecZ;
228 return Input_RecZ+ShiftValue;
229}
230
231void ReadCosmicRayData::ReadCgemDigits()
232{
233 //cout<<"Start to read the CgemDigits"<<endl;
234 Tdigi->GetEntry(Ind_Entry_D);
235 Ind_Entry_D++;
236}
237
238void ReadCosmicRayData::ReadCgemClusters()
239{
240 //cout<<"Start to read the CgemClusters"<<endl;
241 Tcluster->GetEntry(Ind_Entry_C);
242 Ind_Entry_C++;
243}
244
245bool ReadCosmicRayData::ConvertHitToDigi(int ihit, unsigned int &charge_channel, unsigned int &time_channel)
246{
247
248 if(CosmicRayDataSetID == "CR201909") {
249 // cout << "ReadCosmicRayData::ConvertGRAALToCgemBoss(), converting set " << CosmicRayDataSetID << endl;
250
251 // if(!m_GemHit_is_tpc[i]&&Cut_on_tpc) return false; // ??
252 charge_channel = 1;
253 time_channel = 1;
254
255 // from plane to LayerID
256 if(m_plane[ihit] == 0) m_LayerID[ihit] = 0; // LAYER1
257 else if(m_plane[ihit] == 1) m_LayerID[ihit] = 1; // LAYER2
258 else if(m_plane[ihit] == 2) m_LayerID[ihit] = 1; // LAYER2
259 else return false;
260
261 // from view to StripType
262 if(m_view[ihit] == 2) m_StripType[ihit] = 0; // x strip (phi)
263 else if(m_view[ihit] == 3) m_StripType[ihit] = 1; // v strip (stereo)
264 else return false;
265
266 // sheet ID/strip ID
267 m_SheetID[ihit] = 0; // LAYER 1
268 m_StripID[ihit] = m_strip[ihit];
269 if(m_LayerID[ihit] == 1) { // LAYER 2
270 if(m_StripType[ihit] == 0 && m_strip[ihit] >= CgemID::getXSTRIP_MAX(m_LayerID[ihit])) {
271 m_StripID[ihit] = m_strip[ihit] - CgemID::getXSTRIP_MAX(m_LayerID[ihit]);
272 m_SheetID[ihit] = 1;
273 }
274 else if(m_StripType[ihit] == 1 && m_strip[ihit] >= CgemID::getVSTRIP_MAX(m_LayerID[ihit])) {
275 m_StripID[ihit] = m_strip[ihit] - CgemID::getVSTRIP_MAX(m_LayerID[ihit]);
276 m_SheetID[ihit] = 1;
277 }
278 }
279
280 // cout << "GRAAL " << m_strip[ihit] << " view " << m_view[ihit] << " plane " << m_plane[ihit] << endl;
281 // cout << "CGEMB " << m_StripID[ihit] << " view " << m_StripType[ihit] << " plane " << m_LayerID[ihit] << " sheet " << m_SheetID[ihit] << endl;
282
283
284 return true;
285 }
286
287 cout << "ERROR : ReadCosmicRayData::ConvertGRAALToCgemBoss(), the data set " << CosmicRayDataSetID << " is unknown! " << endl;
288 return false;
289
290}
291
292void ReadCosmicRayData::SaveCgemDigits()
293{
294
295 bool printFlag=false;
296 //cout<<"Start to save CgemDigi"<<endl;
297 // cgem digis collection defined in BOSS
298 CgemDigiCol* aCgemDigiCol = new CgemDigiCol;
299
300 //cout<<"nDigi : "<<m_nGemHit<<endl;
301 if (m_nGemHit > 0)
302 {
303 // push back cgem digits to CgemDigiCol in BOSS
304 for(int i=0;i<m_nGemHit;i++)
305 {
306 unsigned int charge_channel;
307 unsigned int time_channel;
308 bool is_converted = ConvertHitToDigi(i, charge_channel, time_channel);
309 if(!is_converted) continue;
310 const Identifier ident = CgemID::strip_id(
311 TranslateDigitLayerID(m_LayerID[i]),
312 TranslateDigitSheetID(m_SheetID[i]),
313 TranslateDigitStripType(m_StripType[i]),
314 TranslateDigitStripID(m_StripID[i],TranslateDigitStripType(m_StripType[i])));
315
316 /**
317 cout << "graal strip " << m_strip[i]
318 << " graal view " << m_view[i]
319 << " stripid " << m_StripID[i]
320 << " striptype " << m_StripType[i]
321 << " strip " << TranslateDigitStripID(m_StripID[i],TranslateDigitStripType(m_StripType[i]))
322 << " identifier " << CgemID::getIntID(TranslateDigitLayerID(m_LayerID[i]),
323 TranslateDigitSheetID(m_SheetID[i]),
324 TranslateDigitStripType(m_StripType[i]),
325 TranslateDigitStripID(m_StripID[i],TranslateDigitStripType(m_StripType[i]))) << endl;;
326 **/
327
328 int layerid = CgemID::layer(ident);
329 int stripid = CgemID::strip(ident);
330 int sheetid = CgemID::sheet(ident);
331 bool is_xstrip = CgemID::is_xstrip(ident);
332
333 if(printFlag)
334 {
335 cout
336 << " digiID=" << ident.get_value()
337 << " layerID=" << CgemID::layer(ident)
338 << " sheetID=" << CgemID::sheet(ident)
339 << " isXstrip="<< CgemID::is_xstrip(ident)
340 << " stripID=" << CgemID::strip(ident)
341 << " time channel=" << time_channel
342 << " charge channel=" << charge_channel << endl;
343 }
344
345 CgemDigi* aCgemDigi = new CgemDigi(ident, time_channel, charge_channel);
346 aCgemDigi->setTime_ns(m_time[i]);
347 aCgemDigi->setCharge_fc(m_charge[i]);
348 aCgemDigiCol->push_back(aCgemDigi);
349
350 } // End of 'for(int i=0;i<nDigi;i++)'
351 } // End of 'if (nDigi > 0)'
352
353 // register CGEM digits collection to TDS
354 StatusCode scCgem = m_evtSvc->registerObject("/Event/Digi/CgemDigiCol", aCgemDigiCol);
355 if(scCgem!=StatusCode::SUCCESS)
356 {
357 cout << "ERROR : ReadCosmicRayData::SaveCgemDigits(), Could not register CGEM digi collection! " << endl;
358 }
359
360 //retrieve CGEM digits from TDS
361 /**
362 SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
363 if(!aDigiCol)
364 cout<<"Could not retrieve CGEM digi collection"<<endl;
365
366 CgemDigiCol::iterator iDigiCol;
367 for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
368 {
369 const Identifier ident = (*iDigiCol)->identify();
370 cout<<" layer: "<<CgemID::layer(ident)
371 <<" sheet: "<<CgemID::sheet(ident)
372 <<" strip: "<<CgemID::strip(ident)
373 <<" charge: "<<(*iDigiCol)->getCharge_fc()
374 <<" time: "<<(*iDigiCol)->getTime_ns()<<endl;
375 }
376 cout<<"end of retrieve CGEM digi collection"<<endl;
377 **/
378}
379
380
381
382void ReadCosmicRayData::SaveCgemClusters()
383{
384
385 bool printFlag=false;
386 //cout<<"Start to save CgemCluster"<<endl;
387
388 // cgem digis collection defined in BOSS
389 RecCgemClusterCol* aRecCgemClusterCol = new RecCgemClusterCol;
390
391 int nCluster = m_nGemCluster; // the number of Clusters in one event
392 if (nCluster > 0)
393 {
394 // push back cgem clusters to RecCgemClusterCol in BOSS
395 for(int i=0;i<nCluster;i++)
396 {
397 RecCgemCluster* aRecCgemCluster = new RecCgemCluster();
398
399 aRecCgemCluster->setclusterid(i);
400 aRecCgemCluster->setlayerid(TranslateClusterLayerID(m_ClusterLayerID[i]));
401 aRecCgemCluster->setsheetid(ClusterSheetID);
402 aRecCgemCluster->setflag(TranslateClusterFlag(m_Flag[i]));
403 aRecCgemCluster->setenergydeposit(m_EnergyDeposit[i]);
404 if(TranslateClusterFlag(m_Flag[i])==0)
405 {
406 aRecCgemCluster->setrecphi(m_Cluster_x[i]);
407 aRecCgemCluster->setrecphi_CC(m_Cluster_x_cc[i]);
408 aRecCgemCluster->setrecphi_mTPC(m_Cluster_x_tpc[i]);
409 }
410 if(TranslateClusterFlag(m_Flag[i])==1)
411 {
412 aRecCgemCluster->setrecv(m_Cluster_x[i]);
413 aRecCgemCluster->setrecv_CC(m_Cluster_x_cc[i]);
414 aRecCgemCluster->setrecv_mTPC(m_Cluster_x_tpc[i]);
415 }
416 aRecCgemCluster->setRecZ(ClusterRecZ); //currently, no Z information is available
417 aRecCgemCluster->setRecZ_CC(ClusterRecZ);
418 aRecCgemCluster->setRecZ_mTPC(ClusterRecZ);
419 //aRecCgemCluster->setRecZ(m_Cluster_z[i]);
420 //aRecCgemCluster->setRecZ_CC(m_Cluster_z_cc[i]);
421 //aRecCgemCluster->setRecZ_mTPC(m_Cluster_z_tpc[i]);
422 aRecCgemCluster->setclusterflag(m_ClusterHitIndex[i][0],m_ClusterHitIndex[i][m_ClusternHit[i]-1]);
423
424 if(printFlag)
425 {
426 cout
427 << " clusterID=" << aRecCgemCluster->getclusterid()
428 << " clusterlayerID=" << aRecCgemCluster->getlayerid()
429 << " clustersheetID=" << aRecCgemCluster->getsheetid()
430 << " flag=" << aRecCgemCluster->getflag()
431 << " energydeposit=" << aRecCgemCluster->getenergydeposit()
432 << " recphi=" << aRecCgemCluster->getrecphi()
433 << " recv=" << aRecCgemCluster->getrecv()
434 << " recZ=" << aRecCgemCluster->getRecZ()
435 << " clusterflagb=" << aRecCgemCluster->getclusterflagb()
436 << " clusterflage=" << aRecCgemCluster->getclusterflage()<< endl;
437 }
438
439 aRecCgemClusterCol->push_back(aRecCgemCluster);
440
441 }
442 }
443
444
445 // register CGEM clusters collection to TDS
446 StatusCode scCgem = m_evtSvc->registerObject("/Event/Recon/RecCgemClusterCol", aRecCgemClusterCol); //Where should I put the vector?
447 if(scCgem!=StatusCode::SUCCESS)
448 {
449 cout << "ERROR : ReadCosmicRayData:::SaveCgemClusters(), Could not register CGEM cluster collection! " << endl;
450 }
451
452}
453
454
456
457 MsgStream log(msgSvc(), name());
458 if(ReadDigi&&!ReadCluster) log << MSG::INFO << "ReadCosmicRayData execute(): "<<Ind_Entry_D+1<<"/"<<No_Entries_D<<" events are finished !" << endreq;
459 if(!ReadDigi&&ReadCluster) log << MSG::INFO << "ReadCosmicRayData execute(): "<<Ind_Entry_C+1<<"/"<<No_Entries_C<<" events are finished !" << endreq;
460 if(ReadDigi&&ReadCluster) log << MSG::INFO << "ReadCosmicRayData execute(): "<<Ind_Entry_C+1<<"/"<<No_Entries_C<<" events are finished !" << endreq;
461
462 //interface to event data service
463 ISvcLocator* svcLocator = Gaudi::svcLocator();
464 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
465 if (sc.isFailure())
466 cout<<"Could not accesss EventDataSvc!"<<endl;
467
468 if(ReadDigi)
469 {
470 DigiEvent* aDigiEvent = new DigiEvent;
471 sc = m_evtSvc->registerObject("/Event/Digi",aDigiEvent);
472 if(sc!=StatusCode::SUCCESS) {
473 cout<< "Could not register DigiEvent" <<endl;
474 }
475
476 ReadCgemDigits();
477 SaveCgemDigits();
478 //cout<<"Ind_Entry_D "<<Ind_Entry_D<<", Max_Ind_Entry_D "<<No_Entries_D<<endl;
479
480 if(Ind_Entry_D==No_Entries_D)
481 {
482 log << MSG::INFO << "scheduling a event processing stop...." << endreq;
483 SmartIF<IEventProcessor> ep(serviceLocator());
484 if (ep) ep->stopRun();
485 }
486
487 }
488 if(ReadCluster)
489 {
490 ReconEvent* aReconEvent = new ReconEvent;
491 sc = m_evtSvc->registerObject("/Event/Recon",aReconEvent);
492 if(sc!=StatusCode::SUCCESS) {
493 cout<< "Could not register ReconEvent" <<endl;
494 }
495
496 ReadCgemClusters();
497 SaveCgemClusters();
498 //cout<<"Ind_Entry_C "<<Ind_Entry_C<<", Max_Ind_Entry_C "<<No_Entries_C<<endl;
499 if(Ind_Entry_C==No_Entries_C)
500 {
501 log << MSG::INFO << "scheduling a event processing stop...." << endreq;
502 SmartIF<IEventProcessor> ep(serviceLocator());
503 if (ep) ep->stopRun();
504 }
505 }
506 return StatusCode::SUCCESS;
507}
508
510 MsgStream log(msgSvc(),name());
511 log << MSG::INFO << "ReadCosmicRayData finalize()" << endreq;
512
513 return StatusCode::SUCCESS;
514}
515
516
517
ObjectVector< RecCgemCluster > RecCgemClusterCol
static int strip(const Identifier &id)
static int sheet(const Identifier &id)
static value_type getXSTRIP_MAX(unsigned int f_layer)
static int layer(const Identifier &id)
static value_type getVSTRIP_MAX(unsigned int f_layer)
static bool is_xstrip(const Identifier &id)
static Identifier strip_id(int f_layer, int f_sheet, int f_strip_type, int f_strip)
ReadCosmicRayData(const std::string &name, ISvcLocator *pSvcLocator)