CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemClusterCreate.cxx
Go to the documentation of this file.
1//Head File//
2#include "CgemClusterCreate/CgemClusterCreate.h"
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/AlgFactory.h"
5#include "GaudiKernel/ISvcLocator.h"
6#include "GaudiKernel/DataSvc.h"
7#include "GaudiKernel/SmartDataPtr.h"
8#include "GaudiKernel/IDataProviderSvc.h"
9#include "GaudiKernel/IDataManagerSvc.h"
10#include "GaudiKernel/PropertyMgr.h"
11#include "GaudiKernel/Bootstrap.h"
12#include "GaudiKernel/IService.h"
13#include "GaudiKernel/INTupleSvc.h"
14#include "GaudiKernel/NTuple.h"
15#include "GaudiKernel/RndmGenerators.h"
16#include "EventModel/EventHeader.h"
17#include "KalFitAlg/helix/Helix.h"
18
19#include "RawEvent/RawDataUtil.h"
20#include "RawDataProviderSvc/IRawDataProviderSvc.h"
21#include "Identifier/CgemID.h"
22#include "Identifier/Identifier.h"
23#include "CgemRawEvent/CgemDigi.h"
24
25
26#include "EvTimeEvent/RecEsTime.h"
27
28#include "McTruth/McParticle.h"
29#include "McTruth/CgemMcHit.h"
30
31#include "BesTimerSvc/IBesTimerSvc.h"
32#include "BesTimerSvc/BesTimerSvc.h"
33#include "BesTimerSvc/BesTimer.h"
34
35#include "ReconEvent/ReconEvent.h"
36#include "CgemRecEvent/RecCgemCluster.h"
37
38#include "CLHEP/Vector/ThreeVector.h"
39
40#include <vector>
41#include <iostream>
42#include <cmath>
43#include "TGraphErrors.h"
44#include "TF1.h"
45
46#define W_pitch 0.66
47
48using namespace std;
49using namespace Event;
50//////////////////////////////////////////////////////////////
51//////////////////////////////////////////////////////////////
52const double a_stero[3] = {(45.94*3.14/180),(31.10*3.14/180),(32.99*3.14/180)};
53const double w_sheet[3] = {549.78,417.097,550.614};
54const double dw_sheet[3] = {549.78,416.26,549.78};
55const double r_layer[3] = {87.51,132.7661,175.2661};
56const double dr_layer[3] = {78.338,123.604,166.104};
57const int l_layer[3] = {532,690,847};
58const int n_sheet[3] = {1,2,2};
59/* all the above parameters were only used in method0 which is at very early stage of CGEM software
60 2016-12-29 wll
61 */
62
63// --- for mTPC fit1
64//const double time_mid = -40.0;
65// -----------------
66
67// --- for mTPC fit2
68//const double time0 = -110;
69//const double drift_velocity= 38/1000.;
70const double drift_gap = 5;
71// -----------------
72
73CgemClusterCreate::CgemClusterCreate(const std::string& name, ISvcLocator* pSvcLocator):
74 Algorithm(name,pSvcLocator)
75{
76 declareProperty("PrintFlag",myPrintFlag=0);
77 declareProperty("MotherParticleID",myMotherParticleID=443);
78 declareProperty("Method",myMethod=2);
79 declareProperty("ntuple",myNtuple=0);
80 declareProperty("effCluster",myClusterEff=1.0);
81 declareProperty("fillOpt",m_fillOption=-1);
82 declareProperty("selGoodDigi",m_selGoodDigi=1);
83 declareProperty("minDigiTime",m_minDigiTime=-8875);
84 declareProperty("maxDigiTime",m_maxDigiTime=-8562);
85 declareProperty("TPCFitMethod",m_selectTPC=1);
86 declareProperty("minQDigi",myQMin=0.0);
87 declareProperty("LUTfile", m_LUTfile = "LUTfile.root");
88 declareProperty("clustering_gap", clustering_gap = 0);
89 declareProperty("dropHitRandomly", myDropHit = 0);
90 declareProperty("uTPC_correction", uTPC_correction = true);
91
92 // --- for random noise (toy cluster)
93 declareProperty("addRandomToyCluster", myRandomCluster=0);
94 declareProperty("NoiseRate_layer1", myNoiseRate[0]=0.0);
95 declareProperty("NoiseRate_layer2", myNoiseRate[1]=0.0);
96 declareProperty("NoiseRate_layer3", myNoiseRate[2]=0.0);
97
98
99 if(myMethod==0) reset();
100 else if(myMethod==1||myMethod==3) resetFiredStripMap();
101
102 m_totEvt = 0;
103 m_nStripRead=0;
104 m_iStripToDrop=-1;
105}
106
107CgemClusterCreate::~CgemClusterCreate()
108{
109 delete lutreader;
110}
111
112StatusCode CgemClusterCreate::initialize()
113{
114 MsgStream log(msgSvc(), name());
115 log << MSG::INFO << "in initialize()" << endreq;
116
117 //CgemGeomSvc//
118 // StatusCode sc;
119 // ICgemGeomSvc* ISvcCgem;
120 // sc = service("CgemGeomSvc", ISvcCgem);
121 // if(sc != StatusCode::SUCCESS)
122 // {
123 // log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
124 // return StatusCode::FAILURE;
125 // }
126
127 //BesTimerSvc//
128 StatusCode tsc;
129 IBesTimerSvc* m_timersvc;
130 tsc = service( "BesTimerSvc", m_timersvc);
131 if( tsc.isFailure() )
132 {
133 log << MSG::WARNING << name() << " Unable to locate BesTimerSvc" << endreq;
134 return StatusCode::FAILURE;
135 }
136 m_timer = m_timersvc->addItem("Execution");
137
138 if(myNtuple) myNTupleHelper=new NTupleHelper(ntupleSvc()->book("RecCgem/CgemCluster",CLID_ColumnWiseTuple,"CgemCluster"));
139
140 if(myMethod==0||myMethod==2) {
141 if(myNtuple) hist_def();
142 }
143
144 // --- get CgemGeomSvc ---
145 ISvcLocator* svcLocator = Gaudi::svcLocator();
146 ICgemGeomSvc* ISvc;
147 StatusCode sc=svcLocator->service("CgemGeomSvc", ISvc);
148 myCgemGeomSvc=dynamic_cast<CgemGeomSvc *>(ISvc);
149 //StatusCode sc = service( "CgemGeomSvc", myCgemGeomSvc);
150 if (!sc.isSuccess()) log<< MSG::INFO << "CgemClusterCreate::initialize(): Could not open CGEM geometry file" << endreq;
151 myNCgemLayers = myCgemGeomSvc->getNumberOfCgemLayer();
152 if(myPrintFlag) cout<<"CgemClusterCreate::initialize() "<<myNCgemLayers<<" Cgem layers"<<endl;
153 for(int i=0; i<myNCgemLayers; i++)
154 {
155 myCgemLayer[i]=myCgemGeomSvc->getCgemLayer(i);
156 myNSheets[i]=myCgemLayer[i]->getNumberOfSheet();
157 bool IsReverse = myCgemLayer[i]->getOrientation();
158 myRXstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu2();
159 myRVstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu1();
160 if(IsReverse) {
161 myRXstrips[i] = myCgemLayer[i]->getOuterROfAnodeCu2();
162 myRVstrips[i] = myCgemLayer[i]->getOuterROfAnodeCu1();
163 }
164 if(myRandomCluster>0) {
165 myNoiseRate[i]*=myCgemLayer[i]->getNumberOfChannelX();
166 cout<<"going to add toy cluster as noise with noise rate "<<myNoiseRate[i]<<" for CGEM layer "<<i<<endl;
167 }
168 if(myPrintFlag) cout<<"layer "<<i<<": "<<myNSheets[i]<<" sheets"<<", is reversed "<<IsReverse<<", RX="<<myRXstrips[i]<<", RY="<<myRVstrips[i]<<endl;
169 }
170
171 sc = service("CgemGeomSvc", m_SvcCgem);
172 if(sc != StatusCode::SUCCESS) {
173 log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
174 return StatusCode::FAILURE;
175 }
176
177 anode = m_SvcCgem->getReadoutPlane(0, 0);
178 anode_mid_gap_radius[0] = anode->getMidRAtGap();
179 anode = m_SvcCgem->getReadoutPlane(1, 0);
180 anode_mid_gap_radius[1] = anode->getMidRAtGap();
181 anode = m_SvcCgem->getReadoutPlane(2, 0);
182 anode_mid_gap_radius[2] = anode->getMidRAtGap();
183
184
185 // --- get CgemCalibFunSvc ---
186 sc = service ("CgemCalibFunSvc", myCgemCalibSvc);
187 if ( sc.isFailure() ){
188 cout<< name() << "Could not load MdcCalibFunSvc!" << endl;
189 return sc;
190 }
191 //cout<<"sigma from CgemCalibFunSvc: "<<myCgemCalibSvc->getSigma(0,0,0,0,0,0)<<endl;
192
193
194 // LUT
195 lutreader = new CgemLUTReader(m_LUTfile);
196 lutreader->ReadLUT();
197
198 return StatusCode::SUCCESS;
199}
200
201StatusCode CgemClusterCreate::execute()
202{
203 MsgStream log(msgSvc(), name());
204 log << MSG::INFO << "in execute()" << endreq;
205
206 // check and register /Event/Recon
207 DataObject *aReconEvent;
208 eventSvc()->findObject("/Event/Recon",aReconEvent);
209 if(aReconEvent==NULL) {
210 aReconEvent = new ReconEvent();
211 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
212 if(sc!=StatusCode::SUCCESS) {
213 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
214 return( StatusCode::FAILURE);
215 }
216 }
217
218
219 if(myMethod==0) method0();
220 else if(myMethod==1) method1();
221 else if(myMethod==2) toyCluster();
222 else if(myMethod==3) method2();
223 m_totEvt++;
224
225 return StatusCode::SUCCESS;
226}
227
228StatusCode CgemClusterCreate::method0()
229{
230 //timer begin//
231 m_timer->start();
232
233 MsgStream log(msgSvc(), name());
234 log << MSG::INFO << "in method0()" << endreq;
235
236 //clear the container//
237 reset();
238
239 // get new information//
240 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
241 if (!evtHead)
242 {
243 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
244 return StatusCode::FAILURE;
245 }
246 int run=evtHead->runNumber();
247 int evt=evtHead->eventNumber();
248 //if(evt!=9180) return StatusCode::SUCCESS;// to check with a specific event
249
250 if(myPrintFlag) {
251 cout<<endl;
252 cout<<"-----------------------------------------------"<<endl;
253 cout<<"CgemClusterCreate::execute: run "<<run<<", evt "<<evt<<endl;
254 }
255
256
257 SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(),"/Event/Digi/CgemDigiCol");
258 if (!cgemDigiCol)
259 {
260 log << MSG::WARNING << "Could not retrieve Cgem digi list" << endreq;
261 return StatusCode::FAILURE;
262 }
263
264 //loop every digi//
265 CgemDigiCol::iterator iter_digi = cgemDigiCol->begin();
266 int layerid,sheetid,stripid,time_chnnel;
267 double energydeposit;
268 double nXStrips[3]={0,0,0};
269 double nVStrips[3]={0,0,0};
270 for( ; iter_digi != cgemDigiCol->end(); ++iter_digi)
271 {
272 layerid = CgemID::layer((*iter_digi)->identify());
273 sheetid = CgemID::sheet((*iter_digi)->identify());
274 stripid = CgemID::strip((*iter_digi)->identify());
275 energydeposit = (*iter_digi)->getChargeChannel();
276 time_chnnel = (*iter_digi)->getTimeChannel();
277 bool striptype = CgemID::is_xstrip((*iter_digi)->identify());
278 int flagxv;
279 if(striptype == true)
280 {
281 flagxv = 0;
282 nXStrips[layerid]+=1;
283 }
284 else {
285 flagxv = 1;
286 nVStrips[layerid]+=1;
287 }
288 if(myPrintFlag)
289 {
290 if(iter_digi==cgemDigiCol->begin())
291 {
292 cout<<"cgemDigiCol:"<<endl;
293 cout
294 <<setw(10)<<"layer"
295 <<setw(10)<<"sheet"
296 <<setw(10)<<"XV(0/1)"
297 <<setw(10)<<"strip_ID"
298 <<setw(15)<<"E"
299 <<setw(15)<<"T"
300 <<endl;
301 }
302 cout
303 <<setw(10)<<layerid
304 <<setw(10)<<sheetid
305 <<setw(10)<<flagxv
306 <<setw(10)<<stripid
307 <<setw(15)<<setprecision(10)<<energydeposit
308 <<setw(15)<<setprecision(10)<<time_chnnel
309 <<endl;
310 }
311 m_strip[layerid][sheetid][flagxv][stripid] = 1;
312 m_edep[layerid][sheetid][flagxv][stripid] = energydeposit;
313
314 }
315
316 processstrips(0);
317 processstrips(1);
318 transcluster();
319 mixcluster();
320
321 RecCgemClusterCol* clustercol = new RecCgemClusterCol;
322
323 //register RecCgemClusterCol//
324 IDataProviderSvc* evtSvc = NULL;
325 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
326 if (evtSvc) {
327 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
328 } else {
329 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
330 return StatusCode::SUCCESS;
331 }
332
333 StatusCode clustersc;
334 IDataManagerSvc *dataManSvc;
335 dataManSvc= dynamic_cast<IDataManagerSvc*>(evtSvc);
336 DataObject *aClusterCol;
337 evtSvc->findObject("/Event/Recon/RecCgemClusterCol",aClusterCol);
338 if(aClusterCol != NULL) {
339 dataManSvc->clearSubTree("/Event/Recon/RecCgemClusterCol");
340 evtSvc->unregisterObject("/Event/Recon/RecCgemClusterCol");
341 }
342
343 clustersc = evtSvc->registerObject("/Event/Recon/RecCgemClusterCol", clustercol);
344 if( clustersc.isFailure() ) {
345 log << MSG::FATAL << "Could not register RecCgemCluster" << endreq;
346 return StatusCode::SUCCESS;
347 }
348 log << MSG::INFO << "RecCgemClusterCol registered successfully!" <<endreq;
349
350 //push back the cluster from the map//
351 for(m_x_map_it = m_x_map.begin();m_x_map_it!=m_x_map.end();++m_x_map_it){
352 RecCgemCluster* reccluster = (*m_x_map_it).second;
353 clustercol->push_back(reccluster);
354 }
355 for(m_v_map_it = m_v_map.begin();m_v_map_it!=m_v_map.end();++m_v_map_it){
356 RecCgemCluster* reccluster = (*m_v_map_it).second;
357 clustercol->push_back(reccluster);
358 }
359 for(m_xv_map_it = m_xv_map.begin();m_xv_map_it!=m_xv_map.end();++m_xv_map_it){
360 RecCgemCluster* reccluster = (*m_xv_map_it).second;
361 clustercol->push_back(reccluster);
362 }
363
364 //read the information to fill the histogram//
365 SmartDataPtr<RecCgemClusterCol> cgemClusterCol(eventSvc(),"/Event/Recon/RecCgemClusterCol");
366 if (!cgemClusterCol){
367 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
368 return StatusCode::FAILURE;
369 }
370
371 m_evt = evtHead->eventNumber();
372 m_evt1 = evtHead->eventNumber();
373
374 // nt1 //
375 RecCgemClusterCol::iterator iter_cluster = cgemClusterCol->begin();
376 int ii=0;
377 for( ; iter_cluster != cgemClusterCol->end(); ++iter_cluster){
378 if((*iter_cluster)->getflag()==2){
379 if(m_nt1){
380 m_layerid[ii] = (*iter_cluster)->getlayerid();
381 m_sheetid[ii] = (*iter_cluster)->getsheetid();
382 m_clusterid[ii]= (*iter_cluster)->getclusterid();
383 m_flag[ii] = (*iter_cluster)->getflag();
384 m_energy[ii] = (*iter_cluster)->getenergydeposit();
385 m_recphi[ii] = (*iter_cluster)->getrecphi();
386 m_recv[ii] = (*iter_cluster)->getrecv();
387 }
388 ii++;
389 }
390 else{continue;}
391 }
392 m_ncluster = ii;
393 m_nt1->write();
394
395 // read the truth information nt3 and nt6 //
396 m_timer->stop();
397 m_evttime=m_timer->elapsed();
398 m_nt5->write();
399
400 return StatusCode::SUCCESS;
401}
402
403StatusCode CgemClusterCreate::method1()
404{
405 //timer begin//
406 m_timer->start();
407
408 MsgStream log(msgSvc(), name());
409 log << MSG::INFO << "in method1()" << endreq;
410
411 //clear the container//
412 resetFiredStripMap();
413
414 // get new information//
415 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
416 if (!evtHead)
417 {
418 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
419 return StatusCode::FAILURE;
420 }
421 int run=evtHead->runNumber();
422 int evt=evtHead->eventNumber();
423 //if(evt!=996) return StatusCode::SUCCESS;// to check with a specific event
424
425 if(myPrintFlag) {
426 cout<<endl;
427 cout<<"-----------------------------------------------"<<endl;
428 cout<<"CgemClusterCreate::execute: run "<<run<<", evt "<<evt<<endl;
429 }
430
431 if(run<0) fillMCTruth();
432
433 SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(),"/Event/Digi/CgemDigiCol");
434 if (!cgemDigiCol)
435 {
436 log << MSG::WARNING << "Could not retrieve Cgem digi list" << endreq;
437 return StatusCode::FAILURE;
438 }
439
440 //loop every digi//
441 CgemDigiCol::iterator iter_digi = cgemDigiCol->begin();
442 int layerid,sheetid,stripid,time_chnnel;
443 double energydeposit;
444 double Q_fC, T_ns;
445 double nXStrips[3]={0,0,0};
446 double nVStrips[3]={0,0,0};
447 bool printed = false;
448 for( ; iter_digi != cgemDigiCol->end(); ++iter_digi)
449 {
450 //if( m_selGoodDigi== 1 && !isGoodDigi(iter_digi) ) continue;
451 //if( m_selGoodDigi==-1 && isGoodDigi(iter_digi) ) continue;
452 if(!selDigi(iter_digi,m_selGoodDigi)) continue;
453 layerid = CgemID::layer((*iter_digi)->identify());
454 sheetid = CgemID::sheet((*iter_digi)->identify());
455 stripid = CgemID::strip((*iter_digi)->identify());
456 energydeposit = (*iter_digi)->getChargeChannel();
457 time_chnnel = (*iter_digi)->getTimeChannel();
458 Q_fC = (*iter_digi)->getCharge_fc();
459 T_ns = (*iter_digi)->getTime_ns();
460 bool striptype = CgemID::is_xstrip((*iter_digi)->identify());
461 int flagxv;
462 if(striptype == true)
463 {
464 flagxv = 0;
465 nXStrips[layerid]+=1;
466 }
467 else {
468 flagxv = 1;
469 nVStrips[layerid]+=1;
470 }
471 if(myPrintFlag)
472 {
473 //if(iter_digi==cgemDigiCol->begin())
474 if(!printed)
475 {
476 cout<<"cgemDigiCol:"<<endl;
477 cout
478 <<setw(10)<<"layer"
479 <<setw(10)<<"sheet"
480 <<setw(10)<<"XV(0/1)"
481 <<setw(10)<<"strip_ID"
482 <<setw(15)<<"E"
483 <<setw(15)<<"T"
484 <<endl;
485 printed=true;
486 }
487 cout
488 <<setw(10)<<layerid
489 <<setw(10)<<sheetid
490 <<setw(10)<<flagxv
491 <<setw(10)<<stripid
492 //<<setw(15)<<setprecision(10)<<energydeposit
493 //<<setw(15)<<setprecision(10)<<time_chnnel
494 <<setw(15)<<setprecision(10)<<Q_fC
495 <<setw(15)<<setprecision(10)<<T_ns
496 <<endl;
497 }
498 myFiredStripMap[layerid][sheetid][flagxv][stripid]=iter_digi;
499 }
500
501 //register RecCgemClusterCol//
502 SmartDataPtr<RecCgemClusterCol> lastCgemClusterCol(eventSvc(),"/Event/Recon/RecCgemClusterCol");
503 IDataProviderSvc* evtSvc = NULL;
504 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
505 if (evtSvc) {
506 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
507 } else {
508 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
509 return StatusCode::SUCCESS;
510 }
511 if(lastCgemClusterCol) {
512 evtSvc->unregisterObject("/Event/Recon/RecCgemClusterCol");
513 //cout<<"RecCgemClusterCol not found"<<endl;
514 }
515 StatusCode sc;
516 RecCgemClusterCol* aCgemClusterCol = new RecCgemClusterCol;
517 sc = evtSvc->registerObject("/Event/Recon/RecCgemClusterCol", aCgemClusterCol);
518 if(sc.isFailure())
519 {
520 log << MSG::FATAL << "Could not register RecCgemClusterCol" << endreq;
521 return StatusCode::SUCCESS;
522 }
523
524 // start cluster reconstruction
525 if(myPrintFlag) cout<<"--------------------------------------------"<<endl;
526 // --- for CC
527 double sumQ(0.0),sumQX(0.0);
528
529 int nCluster[3][3]={{0,0,0},{0,0,0},{0,0,0}};//[layer][XV] (0:X, 1:V, 2:XV)
530 double posX[3][100],posZ[3][100],posV[3][100],phi_XV[3][100];
531 double QX[3][100],QV[3][100],QX_XV[3][100],QV_XV[3][100];
532 int nStripsX[3][100],nStripsV[3][100];
533
534 vector<int> iStart_cluster[3][2][2];//[layer][sheet][XV]
535 vector<int> iEnd_cluster[3][2][2];
536 vector<double> E_cluster[3][2][2];//[layer][sheet][X or V or XV]
537 vector<int> id_cluster[3][2][3];// index in aCgemClusterCol
538 vector<double> vecPos_cluster[3][2][3];//[layer][sheet][X or V or XV]
539 vector<int> idxCluster[3][2][2]; // cluster index in iStart_cluster etc for XV clusters, [layer][sheet][XV]
540 vector<int> idxBoundaryXVcluster[3][2][2];// keep the index in idxCluster for XV cluster at X boundaries, [layer][sheet][end or start]
541
542 //double posX_driftRegion[100];
543 for(int i=0; i<3; i++)//layer
544 {
545 if(myPrintFlag) cout<<"---- layer "<<i<<" ----"<<endl;
546 for(int j=0; j<myNSheets[i]; j++)// sheet
547 {
548 if(myPrintFlag) cout<<"---- sheet "<<j<<":: "<<endl;
549 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(i,j);
550 for(int k=0; k<2; k++) // XV
551 {
552 if(myPrintFlag) cout<<"---- XV "<<k<<": "<<endl;
553
554 map<int,CgemDigiCol::iterator>::iterator iter=myFiredStripMap[i][j][k].begin();
555 map<int,CgemDigiCol::iterator>::iterator iter_end=myFiredStripMap[i][j][k].end();
556 int i1(-1),i2(-1),idLast(-2);
557 for(; iter!=iter_end; iter++)
558 {
559 if((iter->first-idLast)>1)
560 {
561 if(idLast!=-2)// find one cluster
562 {
563 double posCluster(9999.);
564 if(sumQ>0) {
565 posCluster = sumQX/sumQ;
566 if(myPrintFlag) {
567 if(k==0) cout<<"find X cluster "<<nCluster[i][k]<<" from strip "<<i1<<" to "<<i2<<", position = "<<posCluster<<" rad"<<endl;
568 if(k==1) cout<<"find V cluster "<<nCluster[i][k]<<" from strip "<<i1<<" to "<<i2<<", position = "<<posCluster<<" mm"<<endl;
569 }
570 int clusterId=aCgemClusterCol->size();
571 iStart_cluster[i][j][k].push_back(i1);
572 iEnd_cluster[i][j][k].push_back(i2);
573 vecPos_cluster[i][j][k].push_back(posCluster);
574 E_cluster[i][j][k].push_back(sumQ);
575 id_cluster[i][j][k].push_back(clusterId);
576
577 RecCgemCluster *pt_recCluster = new RecCgemCluster();
578 pt_recCluster->setclusterid(clusterId);
579 pt_recCluster->setlayerid(i);
580 pt_recCluster->setsheetid(j);
581 pt_recCluster->setflag(k);
582 pt_recCluster->setenergydeposit(sumQ);
583 pt_recCluster->setclusterflag(i1,i2);
584 if(k==0) pt_recCluster->setrecphi(posCluster);
585 if(k==1) pt_recCluster->setrecv(posCluster);
586 aCgemClusterCol->push_back(pt_recCluster);
587
588 if(k==0&&nCluster[i][k]<100) // X clusters
589 {
590 posX[i][nCluster[i][k]]=posCluster;
591 QX[i][nCluster[i][k]]=sumQ;
592 nStripsX[i][nCluster[i][k]]=i2-i1+1;
593 // correction to drift region
594 }
595 if(k==1)// V clusters
596 {
597 if(nCluster[i][k]<100) {
598 nStripsV[i][nCluster[i][k]]=i2-i1+1;
599 posV[i][nCluster[i][k]]=posCluster;
600 QV[i][nCluster[i][k]]=sumQ;
601 }
602 int nXCluster=iStart_cluster[i][j][0].size();
603 for(int iX=0; iX<nXCluster; iX++)
604 {
605 double Z_pos = readoutPlane->getZFromPhiV(vecPos_cluster[i][j][0][iX],posCluster);
606 //if(Z_pos!=-9999.0) // XV clusters found
607 if(readoutPlane->OnThePlane(vecPos_cluster[i][j][0][iX],Z_pos)) // XV clusters found
608 {
609 if(myPrintFlag) cout<<"find a XV cluster, Z="<<Z_pos<<endl;
610 int iV=iStart_cluster[i][j][1].size()-1;
611 clusterId=aCgemClusterCol->size();
612 vecPos_cluster[i][j][2].push_back(Z_pos);
613 id_cluster[i][j][2].push_back(clusterId);
614 idxCluster[i][j][0].push_back(iX);
615 idxCluster[i][j][1].push_back(iV);
616
617 RecCgemCluster *pt2_recCluster = new RecCgemCluster();
618 pt2_recCluster->setclusterid(clusterId);
619 pt2_recCluster->setlayerid(i);
620 pt2_recCluster->setsheetid(j);
621 pt2_recCluster->setflag(2);
622 pt2_recCluster->setclusterflag(id_cluster[i][j][0][iX],id_cluster[i][j][1][iV]);
623 pt2_recCluster->setrecphi(vecPos_cluster[i][j][0][iX]);
624 pt2_recCluster->setrecv(vecPos_cluster[i][j][1][iV]);
625 pt2_recCluster->setRecZ(Z_pos);
626 pt2_recCluster->setenergydeposit(sumQ+E_cluster[i][j][0][iX]);
627 aCgemClusterCol->push_back(pt2_recCluster);
628
629 if(nCluster[i][2]<100) {
630 posZ[i][nCluster[i][2]]=Z_pos;
631 phi_XV[i][nCluster[i][2]]=vecPos_cluster[i][j][0][iX];
632 }
633 if( (iEnd_cluster[i][j][0][iX]+1)==readoutPlane->getNXstrips() )
634 idxBoundaryXVcluster[i][j][0].push_back(idxCluster[i][j][1].size()-1);
635 if(iStart_cluster[i][j][0][iX]==0)
636 idxBoundaryXVcluster[i][j][1].push_back(idxCluster[i][j][1].size()-1);
637 nCluster[i][2]++;
638 }
639 }
640 }
641 nCluster[i][k]++;
642 }
643 else {
644 cout<<"sumQ<=0!: sumQX,sumQ="<<sumQX<<","<<sumQ<<endl;
645 }
646 sumQ=0.0;
647 sumQX=0.0;
648 }// find one cluster
649 i1=iter->first;
650 i2=iter->first;
651 }
652 else {
653 i2=iter->first;
654 }
655 idLast=iter->first;
656 if(myPrintFlag) cout<<"fired strip "<<idLast<<endl;
657 //double energy=(*(iter->second))->getChargeChannel();
658 double energy=(*(iter->second))->getCharge_fc();
659 double pos=9999.;
660 if(k==0) pos=readoutPlane->getPhiFromXID(iter->first);
661 if(k==1) pos=readoutPlane->getCentralVFromVID(iter->first);
662 sumQ+=energy;
663 sumQX+=pos*energy;
664 }// loop fired strips
665
666 if(idLast!=-2)// find the last cluster
667 {
668 double posCluster(9999.);
669 if(sumQ>0) {
670 posCluster = sumQX/sumQ;
671 if(myPrintFlag) {
672 if(k==0) cout<<"find X cluster "<<nCluster[i][k]<<" from strip "<<i1<<" to "<<i2<<", position = "<<posCluster<<" rad"<<endl;
673 if(k==1) cout<<"find V cluster "<<nCluster[i][k]<<" from strip "<<i1<<" to "<<i2<<", position = "<<posCluster<<" mm"<<endl;
674 }
675 int clusterId=aCgemClusterCol->size();
676 iStart_cluster[i][j][k].push_back(i1);
677 iEnd_cluster[i][j][k].push_back(i2);
678 vecPos_cluster[i][j][k].push_back(posCluster);
679 E_cluster[i][j][k].push_back(sumQ);
680 id_cluster[i][j][k].push_back(clusterId);
681
682 RecCgemCluster *pt_recCluster = new RecCgemCluster();
683 pt_recCluster->setclusterid(clusterId);
684 pt_recCluster->setlayerid(i);
685 pt_recCluster->setsheetid(j);
686 pt_recCluster->setflag(k);
687 pt_recCluster->setenergydeposit(sumQ);
688 pt_recCluster->setclusterflag(i1,i2);
689 if(k==0) pt_recCluster->setrecphi(posCluster);
690 if(k==1) pt_recCluster->setrecv(posCluster);
691 aCgemClusterCol->push_back(pt_recCluster);
692
693 if(k==0&&nCluster[i][k]<100) {
694 posX[i][nCluster[i][k]]=posCluster;
695 QX[i][nCluster[i][k]]=sumQ;
696 nStripsX[i][nCluster[i][k]]=i2-i1+1;
697 // correction to drift region
698 }
699 if(k==1)// V clusters
700 {
701 if(nCluster[i][k]<100) {
702 nStripsV[i][nCluster[i][k]]=i2-i1+1;
703 posV[i][nCluster[i][k]]=posCluster;
704 QV[i][nCluster[i][k]]=sumQ;
705 }
706 int nXCluster=iStart_cluster[i][j][0].size();
707 for(int iX=0; iX<nXCluster; iX++)
708 {
709 double Z_pos = readoutPlane->getZFromPhiV(vecPos_cluster[i][j][0][iX],posCluster);
710 //if(Z_pos!=-9999.0) // XV clusters found
711 if(readoutPlane->OnThePlane(vecPos_cluster[i][j][0][iX],Z_pos)) // XV clusters found
712 {
713 if(myPrintFlag) cout<<"find a XV cluster, Z="<<Z_pos<<endl;
714 clusterId=aCgemClusterCol->size();
715 vecPos_cluster[i][j][2].push_back(Z_pos);
716 id_cluster[i][j][2].push_back(clusterId);
717 int iV=iStart_cluster[i][j][1].size()-1;
718
719 RecCgemCluster *pt2_recCluster = new RecCgemCluster();
720 pt2_recCluster->setclusterid(clusterId);
721 pt2_recCluster->setlayerid(i);
722 pt2_recCluster->setsheetid(j);
723 pt2_recCluster->setflag(2);
724 pt2_recCluster->setclusterflag(id_cluster[i][j][0][iX],id_cluster[i][j][1][iV]);
725 pt2_recCluster->setrecphi(vecPos_cluster[i][j][0][iX]);
726 pt2_recCluster->setrecv(vecPos_cluster[i][j][1][iV]);
727 pt2_recCluster->setRecZ(Z_pos);
728 pt2_recCluster->setenergydeposit(sumQ+E_cluster[i][j][0][iX]);
729 aCgemClusterCol->push_back(pt2_recCluster);
730
731 if(nCluster[i][2]<100) {
732 posZ[i][nCluster[i][2]]=Z_pos;
733 phi_XV[i][nCluster[i][2]]=vecPos_cluster[i][j][0][iX];
734 }
735 idxCluster[i][j][0].push_back(iX);
736 idxCluster[i][j][1].push_back(iEnd_cluster[i][j][1].size()-1);
737 if( (iEnd_cluster[i][j][0][iX]+1)==readoutPlane->getNXstrips() )
738 idxBoundaryXVcluster[i][j][0].push_back(idxCluster[i][j][1].size()-1);
739 if(iStart_cluster[i][j][0][iX]==0)
740 idxBoundaryXVcluster[i][j][1].push_back(idxCluster[i][j][1].size()-1);
741 nCluster[i][2]++;
742 }
743 }
744 }
745 nCluster[i][k]++;
746 }
747 else cout<<"sumQ<=0!: sumQX,sumQ="<<sumQX<<","<<sumQ<<endl;
748 sumQ=0.0;
749 sumQX=0.0;
750 }// find the last cluster on one sheet
751
752 }// loop XV
753 }// loop sheets
754 if(nCluster[i][0]>100) nCluster[i][0]=100;
755 if(nCluster[i][1]>100) nCluster[i][1]=100;
756 if(nCluster[i][2]>100) nCluster[i][2]=100;
757
758 // to combine clusters near the boundary between sheets
759 for(int j=0; j<myNSheets[i]; j++)
760 {
761 // get the index of the next sheet
762 int j_next=j+1;
763 if(j_next==myNSheets[i]) j_next=0;
764
765 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(i,j);
766 CgemGeoReadoutPlane* nextReadoutPlane = myCgemGeomSvc->getReadoutPlane(i,j_next);
767 double xmin_next = nextReadoutPlane->getXmin();
768
769 int nXV_atXEnd=idxBoundaryXVcluster[i][j][0].size();
770 int nXV_atXStart=idxBoundaryXVcluster[i][j_next][1].size();
771 if(nXV_atXEnd==0||nXV_atXStart==0) continue;
772 for(int iXV_atXEnd=0; iXV_atXEnd<nXV_atXEnd; iXV_atXEnd++)
773 {
774 int iXVCluster = idxBoundaryXVcluster[i][j][0][iXV_atXEnd];
775 int iXCluster = idxCluster[i][j][0][iXVCluster];
776 int iVCluster = idxCluster[i][j][1][iXVCluster];
777 int VID1=iStart_cluster[i][j][1][iVCluster];
778 int VID2=iEnd_cluster[i][j][1][iVCluster];
779 int VID1_extend=readoutPlane->getVIDInNextSheetFromVID(VID1, xmin_next);
780 int VID2_extend=readoutPlane->getVIDInNextSheetFromVID(VID2, xmin_next);
781 for(int iXV_atXStart=0; iXV_atXStart<nXV_atXStart; iXV_atXStart++)
782 {
783 int iXVCluster_next = idxBoundaryXVcluster[i][j_next][1][iXV_atXStart];
784 int iXCluster_next = idxCluster[i][j_next][0][iXVCluster_next];
785 int iVCluster_next = idxCluster[i][j_next][1][iXVCluster_next];
786 int VID1_next=iStart_cluster[i][j_next][1][iVCluster_next];
787 int VID2_next=iEnd_cluster[i][j_next][1][iVCluster_next];
788 if( !( (VID1_next-VID2_extend)>1 || (VID1_extend-VID2_next)>1 ) )
789 {
790 int XID1=iStart_cluster[i][j][0][iXCluster];
791 int XID2=iEnd_cluster[i][j][0][iXCluster];
792 int XID1_next=iStart_cluster[i][j_next][0][iXCluster_next];
793 int XID2_next=iEnd_cluster[i][j_next][0][iXCluster_next];
794 int clusterId=aCgemClusterCol->size();
795
796 RecCgemCluster *pt_recCluster = new RecCgemCluster();
797 pt_recCluster->setclusterid(clusterId);
798 pt_recCluster->setlayerid(i);
799 pt_recCluster->setsheetid(-1);
800 pt_recCluster->setflag(3);
801 int id1=id_cluster[i][j][2][iXVCluster];
802 int id2=id_cluster[i][j_next][2][iXVCluster_next];
803 pt_recCluster->setclusterflag(id1,id2);
804 // averaging phi
805 double phi1=vecPos_cluster[i][j][0][iXCluster];
806 double phi2=vecPos_cluster[i][j_next][0][iXCluster_next];
807 if(fabs(phi1-phi2)>CLHEP::pi) // phi domain 0~2pi, so only for 0/2pi
808 {
809 if(phi1>CLHEP::pi) phi1=phi1-CLHEP::twopi;
810 if(phi2>CLHEP::pi) phi2=phi2-CLHEP::twopi;
811 }
812 double E1=E_cluster[i][j][0][iXCluster];
813 double E2=E_cluster[i][j_next][0][iXCluster_next];
814 double phiAverage=(E1*phi1+E2*phi2)/(E1+E2);
815 pt_recCluster->setrecphi(phiAverage);
816 // averging V
817 double v1=vecPos_cluster[i][j][1][iVCluster];
818 v1=readoutPlane->getVInNextSheetFromV(v1, xmin_next);
819 double v2=vecPos_cluster[i][j_next][1][iVCluster_next];
820 E1=E_cluster[i][j][1][iVCluster];
821 E2=E_cluster[i][j_next][1][iVCluster_next];
822 double V_average_next=(E1*v1+E2*v2)/(E1+E2);
823 pt_recCluster->setrecv(V_average_next);
824 // get Z
825 double z_average = nextReadoutPlane->getZFromPhiV(phiAverage,V_average_next,0);
826 pt_recCluster->setRecZ(z_average);
827 // push back into CgemClusterCol
828 aCgemClusterCol->push_back(pt_recCluster);
829
830 if(myPrintFlag) {
831 cout<<"one combinational cluster found: "<<endl;
832 cout<<" sheet "<<j <<" xID:"<<XID1 <<"~"<<XID2 <<", phi:"<<vecPos_cluster[i][j][0][iXCluster] <<", vID:"<<VID1 <<"~"<<VID2 <<", v:"<<vecPos_cluster[i][j][1][iVCluster] <<", z:"<<vecPos_cluster[i][j][2][iXVCluster] <<endl;
833 cout<<" sheet "<<j_next<<" xID:"<<XID1_next<<"~"<<XID2_next<<", phi:"<<vecPos_cluster[i][j_next][0][iXCluster_next]<<", vID:"<<VID1_next<<"~"<<VID2_next<<", v:"<<vecPos_cluster[i][j_next][1][iVCluster_next]<<", z:"<<vecPos_cluster[i][j_next][2][iXVCluster_next]<<endl;
834 cout<<" averagePhi:"<<phiAverage<<", v_average:"<<V_average_next<<", z_average:"<<z_average<<endl;
835 }
836 }
837 }
838 }
839 }// combination of XV clusters
840
841
842 }// loop layers
843
844 // count time used
845 m_timer->stop();
846 double evttime=m_timer->elapsed();
847
848 // checkRecCgemClusterCol
849 if(myPrintFlag) checkRecCgemClusterCol();
850
851 if(myNtuple)
852 {
853 myNTupleHelper->fillLong("run",(long)run);
854 myNTupleHelper->fillLong("evt",(long)evt);
855
856 myNTupleHelper->fillArray("nXstrips","nLayer",(double*) nXStrips,3);
857 myNTupleHelper->fillArray("nVstrips","nLayer",(double*) nVStrips,3);
858
859 myNTupleHelper->fillArray("phiClusterLay1","nXClusterLay1",(double*) posX[0],nCluster[0][0]);
860 myNTupleHelper->fillArrayInt("nXStripsLay1","nXClusterLay1",(int*) nStripsX[0],nCluster[0][0]);
861 myNTupleHelper->fillArray("QXLay1","nXClusterLay1",(double*) QX[0],nCluster[0][0]);
862
863 myNTupleHelper->fillArray("VClusterLay1","nVClusterLay1",(double*) posV[0],nCluster[0][1]);
864 myNTupleHelper->fillArrayInt("nVStripsLay1","nVClusterLay1",(int*) nStripsV[0],nCluster[0][1]);
865 myNTupleHelper->fillArray("QVLay1","nVClusterLay1",(double*) QV[0],nCluster[0][1]);
866
867 myNTupleHelper->fillArray("zClusterLay1","nXVClusterLay1",(double*) posZ[0],nCluster[0][2]);
868 myNTupleHelper->fillArray("phiXVClusterLay1","nXVClusterLay1",(double*) phi_XV[0],nCluster[0][2]);
869
870 myNTupleHelper->fillArray("phiClusterLay2","nXClusterLay2",(double*) posX[1],nCluster[1][0]);
871 myNTupleHelper->fillArrayInt("nXStripsLay2","nXClusterLay2",(int*) nStripsX[1],nCluster[1][0]);
872 myNTupleHelper->fillArray("QXLay2","nXClusterLay2",(double*) QX[1],nCluster[1][0]);
873
874 myNTupleHelper->fillArray("VClusterLay2","nVClusterLay2",(double*) posV[1],nCluster[1][1]);
875 myNTupleHelper->fillArrayInt("nVStripsLay2","nVClusterLay2",(int*) nStripsV[1],nCluster[1][1]);
876 myNTupleHelper->fillArray("QVLay2","nVClusterLay2",(double*) QV[1],nCluster[1][1]);
877
878 myNTupleHelper->fillArray("zClusterLay2","nXVClusterLay2",(double*) posZ[1],nCluster[1][2]);
879 myNTupleHelper->fillArray("phiXVClusterLay2","nXVClusterLay2",(double*) phi_XV[1],nCluster[1][2]);
880
881 myNTupleHelper->fillArray("phiClusterLay3","nXClusterLay3",(double*) posX[2],nCluster[2][0]);
882 myNTupleHelper->fillArrayInt("nXStripsLay3","nXClusterLay3",(int*) nStripsX[2],nCluster[2][0]);
883 myNTupleHelper->fillArray("QXLay3","nXClusterLay3",(double*) QX[2],nCluster[2][0]);
884
885 myNTupleHelper->fillArray("VClusterLay3","nVClusterLay3",(double*) posV[2],nCluster[2][1]);
886 myNTupleHelper->fillArrayInt("nVStripsLay3","nVClusterLay3",(int*) nStripsV[2],nCluster[2][1]);
887 myNTupleHelper->fillArray("QVLay3","nVClusterLay3",(double*) QV[2],nCluster[2][1]);
888
889 myNTupleHelper->fillArray("zClusterLay3","nXVClusterLay3",(double*) posZ[2],nCluster[2][2]);
890 myNTupleHelper->fillArray("phiXVClusterLay3","nXVClusterLay3",(double*) phi_XV[2],nCluster[2][2]);
891
892 myNTupleHelper->fillDouble("evtTime",evttime);
893
894 myNTupleHelper->write();
895 }
896 if(myPrintFlag) cout<<"End of CgemClusterCreate::method1()"<<endl;
897
898 return StatusCode::SUCCESS;
899}// method1
900
901StatusCode CgemClusterCreate::method2()
902{
903 //timer begin//
904 m_timer->start();
905
906 MsgStream log(msgSvc(), name());
907 log << MSG::INFO << "in method1()" << endreq;
908
909 //clear the container//
910 resetFiredStripMap();
911 //resetForCluster();
912
913 // get new information//
914 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
915 if (!evtHead)
916 {
917 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
918 return StatusCode::FAILURE;
919 }
920 int run=evtHead->runNumber();
921 int evt=evtHead->eventNumber();
922 //cout<<"evt = "<<evt<<endl;
923 //if(evt!=996) return StatusCode::SUCCESS;// to check with a specific event
924
925 if(myPrintFlag) {
926 cout<<endl;
927 cout<<"-----------------------------------------------"<<endl;
928 cout<<"CgemClusterCreate::execute: run "<<run<<", evt "<<evt<<endl;
929 }
930
931 if(run<0) fillMCTruth();
932
933 SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(),"/Event/Digi/CgemDigiCol");
934 if (!cgemDigiCol)
935 {
936 log << MSG::WARNING << "Could not retrieve Cgem digi list" << endreq;
937 return StatusCode::FAILURE;
938 }
939
940 //loop every digi//
941 CgemDigiCol::iterator iter_digi = cgemDigiCol->begin();
942 int layerid,sheetid,stripid,time_chnnel;
943 double energydeposit;
944 double Q_fC, T_ns;
945 double nXStrips[3]={0,0,0};
946 double nVStrips[3]={0,0,0};
947 bool printed = false;
948 for( ; iter_digi != cgemDigiCol->end(); ++iter_digi)
949 {
950 if(run<0&&myDropHit==1) {
951 if(m_nStripRead==0) {
952 Rndm::Numbers flat(randSvc(), Rndm::Flat(0,64));
953 m_iStripToDrop=int(flat());
954 //if(myPrintFlag)
955 //cout<<"to drop "<<m_iStripToDrop<<", m_nStripRead=";
956 }
957 //if(myPrintFlag)
958 //cout<<m_nStripRead<<" ";
959 int iStrip = m_nStripRead++;
960 if(m_nStripRead==64) {
961 //if(myPrintFlag)
962 //if(iStrip!=m_iStripToDrop) cout<<endl;
963 m_nStripRead=0;
964 }
965 if(iStrip==m_iStripToDrop) {
966 //if(myPrintFlag)
967 //cout<<"(dropped) ";
968 //if(m_iStripToDrop==63) cout<<endl;
969 continue;
970 }
971 }
972 //if( m_selGoodDigi== 1 && !isGoodDigi(iter_digi) ) continue;
973 //if( m_selGoodDigi==-1 && isGoodDigi(iter_digi) ) continue;
974 if(!selDigi(iter_digi,m_selGoodDigi)) continue;
975 layerid = CgemID::layer((*iter_digi)->identify());
976 sheetid = CgemID::sheet((*iter_digi)->identify());
977 stripid = CgemID::strip((*iter_digi)->identify());
978 energydeposit = (*iter_digi)->getChargeChannel();
979 time_chnnel = (*iter_digi)->getTimeChannel();
980 Q_fC = (*iter_digi)->getCharge_fc();
981 T_ns = (*iter_digi)->getTime_ns();
982 bool striptype = CgemID::is_xstrip((*iter_digi)->identify());
983 int flagxv;
984 if(striptype == true)
985 {
986 flagxv = 0;
987 nXStrips[layerid]+=1;
988 }
989 else {
990 flagxv = 1;
991 nVStrips[layerid]+=1;
992 }
993 if(myPrintFlag)
994 {
995 //if(iter_digi==cgemDigiCol->begin())
996 if(!printed)
997 {
998 cout<<"cgemDigiCol:"<<endl;
999 cout
1000 <<setw(10)<<"layer"
1001 <<setw(10)<<"sheet"
1002 <<setw(10)<<"XV(0/1)"
1003 <<setw(10)<<"strip_ID"
1004 <<setw(15)<<"E"
1005 <<setw(15)<<"T"
1006 <<endl;
1007 printed=true;
1008 }
1009 cout
1010 <<setw(10)<<layerid
1011 <<setw(10)<<sheetid
1012 <<setw(10)<<flagxv
1013 <<setw(10)<<stripid
1014 //<<setw(15)<<setprecision(10)<<energydeposit
1015 //<<setw(15)<<setprecision(10)<<time_chnnel
1016 <<setw(15)<<setprecision(10)<<Q_fC
1017 <<setw(15)<<setprecision(10)<<T_ns
1018 <<endl;
1019 }
1020 myFiredStripMap[layerid][sheetid][flagxv][stripid]=iter_digi;
1021
1022 }
1023
1024 //register RecCgemClusterCol//
1025 SmartDataPtr<RecCgemClusterCol> lastCgemClusterCol(eventSvc(),"/Event/Recon/RecCgemClusterCol");
1026 IDataProviderSvc* evtSvc = NULL;
1027 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
1028 if (evtSvc) {
1029 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1030 } else {
1031 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1032 return StatusCode::SUCCESS;
1033 }
1034 if(lastCgemClusterCol) {
1035 evtSvc->unregisterObject("/Event/Recon/RecCgemClusterCol");
1036 //cout<<"RecCgemClusterCol not found"<<endl;
1037 }
1038 StatusCode sc;
1039 RecCgemClusterCol* aCgemClusterCol = new RecCgemClusterCol;
1040 sc = evtSvc->registerObject("/Event/Recon/RecCgemClusterCol", aCgemClusterCol);
1041 if(sc.isFailure())
1042 {
1043 log << MSG::FATAL << "Could not register RecCgemClusterCol" << endreq;
1044 return StatusCode::SUCCESS;
1045 }
1046
1047 // start cluster reconstruction
1048 if(myPrintFlag) cout<<"--------------------------------------------"<<endl;
1049
1050 // --- for CC
1051 double sumQ(0.0),sumQX(0.0);
1052
1053 // --- for mTPC fit1
1054 vector<double> pos_strips;
1055 vector<double> drift_strips;
1056 vector<double> q_strips;
1057 vector<int> id_strips;
1058 double tposX[3][100],tposZ[3][100],tposV[3][100];//Calculated position by micro-TPC
1059 vector<double> vecTPos_cluster[3][2][3];
1060 // -----------------
1061
1062 // --- for mTPC fit2
1063 double pos_tpc(9999.0), a_tpc(0.0), b_tpc(0.0);
1064 double sum_x_tpc(0.), sum_y_tpc(0.), sum_xx_tpc(0.), sum_yy_tpc(0.), sum_xy_tpc(0.);
1065 int n_tpc(0);
1066 double a_tpc_cluster_X[3][100];
1067 double b_tpc_cluster_X[3][100];
1068 double pos_tpc_cluster_X[3][100];
1069 double a_tpc_cluster_V[3][100];
1070 double b_tpc_cluster_V[3][100];
1071 double pos_tpc_cluster_V[3][100];
1072 // -----------------
1073
1074 int nCluster[3][3]={{0,0,0},{0,0,0},{0,0,0}};//[layer][XV] (0:X, 1:V, 2:XV)
1075 double posX[3][100],posZ[3][100],posV[3][100],phi_XV[3][100];
1076 double QX[3][100],QV[3][100],QX_XV[3][100],QV_XV[3][100];
1077 int nStripsX[3][100],nStripsV[3][100];
1078 vector<int> iStart_cluster[3][2][2];//[layer][sheet][XV]
1079 vector<int> iEnd_cluster[3][2][2];
1080 vector<double> E_cluster[3][2][2];//[layer][sheet][X or V or XV]
1081 vector<int> id_cluster[3][2][3];// index in aCgemClusterCol
1082 vector<double> vecPos_cluster[3][2][3];//[layer][sheet][X or V or XV]
1083 vector<int> idxCluster[3][2][2]; // cluster index in iStart_cluster etc for XV clusters, [layer][sheet][XV]
1084 vector<int> idxBoundaryXVcluster[3][2][2];// keep the index in idxCluster for XV cluster at X boundaries, [layer][sheet][end or start]
1085
1086 //double posX_driftRegion[100];
1087 for(int i=0; i<3; i++)//layer
1088 {
1089 if(myPrintFlag) cout<<"---- layer "<<i<<" ----"<<endl;
1090 for(int j=0; j<myNSheets[i]; j++)// sheet
1091 {
1092 if(myPrintFlag) cout<<"---- sheet "<<j<<":: "<<endl;
1093 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(i,j);
1094 for(int k=0; k<2; k++) // XV
1095 {
1096 if(myPrintFlag) cout<<"---- XV "<<k<<": "<<endl;
1097
1098 map<int,CgemDigiCol::iterator>::iterator iter=myFiredStripMap[i][j][k].begin();
1099 map<int,CgemDigiCol::iterator>::iterator iter_end=myFiredStripMap[i][j][k].end();
1100 map<int,CgemDigiCol::iterator>::iterator iter_next=iter;
1101 int i1(-1),i2(-1),idLast(-2);
1102 bool clusterFound=true;
1103 for(; iter!=iter_end; iter++)
1104 {
1105 if(myPrintFlag) cout<<"fired strip "<<iter->first<<endl;
1106
1107 // get digi charge, time, position
1108 double energy=(*(iter->second))->getCharge_fc();
1109 //double energy=(*(iter->second))->getChargeChannel();
1110 double time=(*(iter->second))->getTime_ns();
1111 double pos=9999.;
1112 if(k==0) pos=readoutPlane->getPhiFromXID(iter->first);
1113 if(k==1) pos=readoutPlane->getCentralVFromVID(iter->first);
1114
1115 // for charge centroid
1116 sumQ+=energy;
1117 sumQX+=pos*energy;
1118
1119 // --- for mTPC fit1
1120 pos_strips.push_back(pos);
1121 // getTimeRising(int layer, int xvFlag, int sheet, int stripID, double Q=100., double z=0.)
1122 double time_rising = myCgemCalibSvc->getTimeRising(i, k, j, iter->first, energy, 0.);
1123 double time_falling = myCgemCalibSvc->getTimeFalling(i, k, j, iter->first, energy, 0.);
1124 double time_gap = time_falling-time_rising;
1125 double drift_velocity = drift_gap/time_gap;
1126 //double time_ns = time-time_rising;
1127 double time_ns=get_Time(iter->second);
1128 if(myPrintFlag) cout<<"T_r, T_f = "<<time_rising<<", "<<time_falling<<", time_ns="<<time<<endl;
1129 drift_strips.push_back(time_ns*drift_velocity);
1130 q_strips.push_back(energy);
1131 id_strips.push_back(iter->first);
1132 // -----------------
1133
1134 // --- for mTPC fit2
1135 double drift = time_ns*drift_velocity;// (time-time0)*drift_velocity;
1136 sum_x_tpc+=pos;
1137 sum_y_tpc+=drift;
1138 sum_xx_tpc+=pos*pos;
1139 sum_yy_tpc+=drift*drift;
1140 sum_xy_tpc+=pos*drift;
1141 n_tpc++;
1142 // -----------------
1143
1144 // check if one cluster found
1145 if(clusterFound) {
1146 i1=iter->first;
1147 clusterFound=false;
1148 }
1149 i2=iter->first;
1150 iter_next++;
1151 int strip_gap_deadstrip = clustering_gap+1;
1152 if(iter_next==iter_end||(iter_next->first-iter->first)>strip_gap_deadstrip) // one cluster found
1153 {
1154 // for charge centroid
1155 double posCluster_CC(9999.);
1156 if(sumQ<=0){
1157 cout<<"sumQ<=0!: sumQX,sumQ="<<sumQX<<","<<sumQ<<endl;
1158 }
1159 else {
1160 posCluster_CC = sumQX/sumQ;
1161 }
1162
1163 // --- for mTPC fit2
1164 double slope2(-9999);
1165 double inter2(-9999);
1166 if(n_tpc>1) {
1167 a_tpc = ((sum_y_tpc*sum_xx_tpc)-(sum_x_tpc*sum_xy_tpc))/(n_tpc*sum_xx_tpc-sum_x_tpc*sum_x_tpc);
1168 b_tpc = (n_tpc*sum_xy_tpc-sum_x_tpc*sum_y_tpc)/(n_tpc*sum_xx_tpc-sum_x_tpc*sum_x_tpc);
1169 if(abs(b_tpc) > 1e-10) { // fix, != 0 is too strict
1170 pos_tpc = (0.5*drift_gap-a_tpc)/b_tpc;
1171 slope2=b_tpc;
1172 inter2=a_tpc;
1173 }
1174 //cout<<"Cluster: "<<aCgemClusterCol->size()<<" "<<n_tpc<<" "<<a_tpc<<" "<<b_tpc<<" "<<pos_tpc<<endl;
1175 }
1176 // ---------------
1177
1178
1179 // --- for mTPC fit1 + output fit1
1180 double tposCluster(9999.);
1181 double slope(-9999);
1182 double inter(-9999);
1183 int n_Strips=pos_strips.size();
1184 if(n_Strips>=2)
1185 {
1186 double *x = new double[n_Strips];
1187 double *tx = new double[n_Strips];
1188 double *ex = new double[n_Strips];
1189 double *etx = new double[n_Strips];
1190 int *stripID = new int[n_Strips];
1191 double x_min = 999;
1192 double x_max = -999;
1193
1194 for(int ix=0; ix<n_Strips; ix++)
1195 {
1196 x[ix] = pos_strips[ix];
1197 tx[ix] = drift_strips[ix];
1198 //ex[ix] = 0.0;
1199 //etx[ix] = 2.0;
1200 ex[ix]=sqrt(pow((W_pitch/anode_mid_gap_radius[i])/sqrt(12.),2)+pow((W_pitch/anode_mid_gap_radius[i])/sqrt(12.)*sumQ/(n_Strips*q_strips[ix]),2));
1201 etx[ix] = 5*drift_velocity;
1202 stripID[ix] = id_strips[ix];
1203 if(x[ix]>x_max) x_max=x[ix];
1204 if(x[ix]<x_min) x_min=x[ix];
1205 if(myPrintFlag)
1206 cout<<"t = "<<tx[ix]<<", q = "<<q_strips[ix]<<", x= "<<x[ix]<<endl;
1207 }
1208
1209 float clsize_real = 1+ stripID[n_Strips-1]-stripID[0];
1210 //Correction
1211 float p0_tpc_corr = 0;
1212 float p1_tpc_corr = 0;
1213 float shift0_tpc_corr = 0;
1214 if(uTPC_correction){
1215 //first round
1216 p0_tpc_corr = 0.003206 + clsize_real * (-0.0006184);
1217 p1_tpc_corr = -0.000500 + clsize_real * ( 0.0001325);
1218 shift0_tpc_corr = 0.003489 + clsize_real * ( 0.0003489);
1219 //second round
1220 p0_tpc_corr += 0.001030 + clsize_real * (-0.0006295);
1221 p1_tpc_corr += 0.000902;
1222 shift0_tpc_corr += 0.001105 + clsize_real * (-0.0002084);
1223 //first rotation
1224 p1_tpc_corr += 0.0005;
1225 //second rotation
1226 p1_tpc_corr += 0.0005;
1227 //third rotation
1228 p1_tpc_corr += 0.0005;
1229 //forth rotation
1230 p1_tpc_corr += 0.0005;
1231 //fifth rotation
1232 p1_tpc_corr += 0.0005;
1233 //sixth rotation
1234 p1_tpc_corr += 0.0010;
1235 //seventh rotation
1236 p1_tpc_corr += 0.0010;
1237 }
1238 vector <double> v_deltaX;
1239 if(k==0){
1240 float real_ix=0;
1241 for(int ix=0; ix<n_Strips; ix++){
1242 float deltaX_corr = 0;
1243 //↗
1244 if(a_tpc*posCluster_CC<0){
1245 if(ix==0){
1246 real_ix=0;
1247 deltaX_corr -= (p0_tpc_corr - shift0_tpc_corr);
1248 }
1249 else{
1250 real_ix+=stripID[ix]-stripID[ix-1];
1251 }
1252 }
1253 //↖
1254 else{
1255 if(ix==n_Strips-1){
1256 real_ix=n_Strips-1;
1257 deltaX_corr += (p0_tpc_corr - shift0_tpc_corr);
1258 }
1259 else{
1260 if(ix!=0) real_ix+=stripID[ix]-stripID[ix-1];
1261 }
1262 }
1263 x[ix] += deltaX_corr;
1264 v_deltaX.push_back(deltaX_corr);
1265 }
1266 }
1267
1268 TGraphErrors xline(n_Strips, x, tx, ex, etx);
1269 //TF1 f1("f1", "[0]*x + [1]", -10000, 10000);
1270 TF1 f1("f1", "[0]*x + [1]", x_min, x_max);
1271 f1.SetParameters(b_tpc,a_tpc); //here
1272 xline.Fit(&f1, "Q");
1273 TF1* f2 = xline.GetFunction("f1");
1274 double xpar[2]={0};
1275 f2->GetParameters(xpar);
1276 double *expar=f2->GetParErrors();
1277 //cout<<xpar[0]<<" +/- "<<expar[0]<<" "<<xpar[1]<<" +/- "<<expar[1]<<endl;
1278 if(xpar[0]!=0)
1279 {
1280 slope=xpar[0];
1281 inter=xpar[1];
1282 tposCluster=(drift_gap/2.-xpar[1])/xpar[0];
1283 if(fabs(tposCluster)>9999) tposCluster=9999.0;
1284 //cout<<"luxl test for mirco tposCluster= "<<tposCluster<<endl;
1285 }
1286 delete[] x;
1287 delete[] tx;
1288 delete[] ex;
1289 delete[] etx;
1290 }
1291 // -----------------
1292 // --- for mTPC fit2
1293 //double slope2(-9999);
1294 //if(n_tpc>1) {
1295 // a_tpc = ((sum_y_tpc*sum_xx_tpc)-(sum_x_tpc*sum_xy_tpc))/(n_tpc*sum_xx_tpc-sum_x_tpc*sum_x_tpc);
1296 // b_tpc = (n_tpc*sum_xy_tpc-sum_x_tpc*sum_y_tpc)/(n_tpc*sum_xx_tpc-sum_x_tpc*sum_x_tpc);
1297 // if(abs(b_tpc) > 1e-10) { // fix, != 0 is too strict
1298 // pos_tpc = (0.5*drift_gap-a_tpc)/b_tpc;
1299 // slope2=b_tpc;
1300 // }
1301 // //cout<<"Cluster: "<<aCgemClusterCol->size()<<" "<<n_tpc<<" "<<a_tpc<<" "<<b_tpc<<" "<<pos_tpc<<endl;
1302 //}
1303 // -----------------
1304
1305 // cout<<"find one cluster!"<<endl;
1306 if(myPrintFlag) {
1307 if(k==0) cout<<"find X cluster "<<nCluster[i][k]<<" from strip "<<i1<<" to "<<i2<<", CC position = "<<posCluster_CC<<" rad"<<", mTPC1 position = "<<tposCluster<<", mTPC2 position = "<<pos_tpc<<endl;
1308 if(k==1) cout<<"find V cluster "<<nCluster[i][k]<<" from strip "<<i1<<" to "<<i2<<", CC position = "<<posCluster_CC<<" mm"<<", mTPC1 position = "<<tposCluster<<", mTPC2 position = "<<pos_tpc<<endl;
1309 //cout<<"find one cluster from strip "<<i1<<" to "<<i2<<", position="<<posCluster_CC<<endl;
1310 }
1311
1312 // fill vectors
1313 int clusterId=aCgemClusterCol->size();
1314 iStart_cluster[i][j][k].push_back(i1);
1315 iEnd_cluster[i][j][k].push_back(i2);
1316 vecPos_cluster[i][j][k].push_back(posCluster_CC);
1317 if(m_selectTPC==1) {
1318 vecTPos_cluster[i][j][k].push_back(tposCluster);
1319 }
1320 else if(m_selectTPC==2) {
1321 vecTPos_cluster[i][j][k].push_back(pos_tpc);
1322 }
1323
1324 E_cluster[i][j][k].push_back(sumQ);
1325 id_cluster[i][j][k].push_back(clusterId);
1326
1327
1328 // --- fill 1D RecCgemCluster ---
1329 RecCgemCluster *pt_recCluster = new RecCgemCluster();
1330 pt_recCluster->setclusterid(clusterId);
1331 pt_recCluster->setlayerid(i);
1332 pt_recCluster->setsheetid(j);
1333 pt_recCluster->setflag(k);
1334 pt_recCluster->setenergydeposit(sumQ);
1335 pt_recCluster->setclusterflag(i1,i2);
1336 if(k==0) {
1337 pt_recCluster->setrecphi(posCluster_CC);
1338 pt_recCluster->setrecphi_CC(posCluster_CC);
1339 if(m_selectTPC==1) {
1340 pt_recCluster->setrecphi_mTPC(tposCluster);
1341 pt_recCluster->setSlope_mTPC(slope);
1342 pt_recCluster->setInter_mTPC(inter);
1343 }
1344 else if(m_selectTPC==2) {
1345 pt_recCluster->setrecphi_mTPC(pos_tpc);
1346 pt_recCluster->setSlope_mTPC(slope2);
1347 pt_recCluster->setInter_mTPC(inter);
1348 }
1349 }
1350 if(k==1) {
1351 pt_recCluster->setrecv(posCluster_CC);
1352 pt_recCluster->setrecv_CC(posCluster_CC);
1353 if(m_selectTPC==1) {
1354 pt_recCluster->setrecv_mTPC(tposCluster);
1355 pt_recCluster->setSlope_mTPC(slope);
1356 pt_recCluster->setInter_mTPC(inter);
1357 }
1358 else if(m_selectTPC==2){
1359 pt_recCluster->setrecv_mTPC(pos_tpc);
1360 pt_recCluster->setSlope_mTPC(slope2);
1361 pt_recCluster->setInter_mTPC(inter2);
1362 }
1363 }
1364 aCgemClusterCol->push_back(pt_recCluster);
1365
1366 // ------------------------------
1367 // --- fill arrays for X-clusters ---
1368 if(k==0&&nCluster[i][k]<100) // if X clusters
1369 {
1370 // --- for CC
1371 nStripsX[i][nCluster[i][k]]=i2-i1+1;
1372 posX[i][nCluster[i][k]]=posCluster_CC;
1373 QX[i][nCluster[i][k]]=sumQ;
1374 // ----------
1375 // --- for mTPC fit1
1376 tposX[i][nCluster[i][k]]=tposCluster;
1377 // -----------------
1378 // --- for mTPC fit2
1379 a_tpc_cluster_X[i][nCluster[i][k]]=a_tpc;
1380 b_tpc_cluster_X[i][nCluster[i][k]]=b_tpc;
1381 pos_tpc_cluster_X[i][nCluster[i][k]]=pos_tpc;
1382 // -----------------
1383 // correction to drift region
1384 }
1385 // -----------------------------------
1386 if(k==1)// if V clusters
1387 {
1388 if(nCluster[i][k]<100) {
1389 // --- for CC
1390 nStripsV[i][nCluster[i][k]]=i2-i1+1;
1391 posV[i][nCluster[i][k]]=posCluster_CC;
1392 QV[i][nCluster[i][k]]=sumQ;
1393 // ----------
1394 // --- for mTPC fit1
1395 tposV[i][nCluster[i][k]]=tposCluster;
1396 // -----------------
1397 }
1398 // --- loop X-V combinations --->
1399 int nXCluster=iStart_cluster[i][j][0].size();
1400 for(int iX=0; iX<nXCluster; iX++)
1401 {
1402 double Z_pos = readoutPlane->getZFromPhiV(vecPos_cluster[i][j][0][iX],posCluster_CC);
1403 double tZ_pos =9999.0;
1404 if(vecTPos_cluster[i][j][0][iX]!=9999.0&&tposCluster!=9999.0)
1405 tZ_pos = readoutPlane->getZFromPhiV(vecTPos_cluster[i][j][0][iX],tposCluster);
1406 //if(Z_pos!=-9999.0) // XV clusters found
1407 if(readoutPlane->OnThePlane(vecPos_cluster[i][j][0][iX],Z_pos)) // XV clusters found
1408 {
1409 if(myPrintFlag) cout<<"find a XV cluster, Z="<<Z_pos<<endl;
1410 int iV=iStart_cluster[i][j][1].size()-1;
1411 clusterId=aCgemClusterCol->size();
1412 vecPos_cluster[i][j][2].push_back(Z_pos);
1413 id_cluster[i][j][2].push_back(clusterId);
1414 idxCluster[i][j][0].push_back(iX);
1415 idxCluster[i][j][1].push_back(iV);
1416 // --- fill 2D RecCgemCluster --->
1417 RecCgemCluster *pt2_recCluster = new RecCgemCluster();
1418 pt2_recCluster->setclusterid(clusterId);
1419 pt2_recCluster->setlayerid(i);
1420 pt2_recCluster->setsheetid(j);
1421 pt2_recCluster->setflag(2);
1422 pt2_recCluster->setenergydeposit(sumQ+E_cluster[i][j][0][iX]);
1423 pt2_recCluster->setclusterflag(id_cluster[i][j][0][iX],id_cluster[i][j][1][iV]);
1424 pt2_recCluster->setrecphi(vecPos_cluster[i][j][0][iX]);
1425 pt2_recCluster->setrecphi_CC(vecPos_cluster[i][j][0][iX]);
1426 pt2_recCluster->setrecphi_mTPC(vecTPos_cluster[i][j][0][iX]);
1427 pt2_recCluster->setrecv(vecPos_cluster[i][j][1][iV]);
1428 pt2_recCluster->setrecv_CC(vecPos_cluster[i][j][1][iV]);
1429 pt2_recCluster->setrecv_mTPC(vecTPos_cluster[i][j][1][iV]);
1430 pt2_recCluster->setRecZ(Z_pos);
1431 pt2_recCluster->setRecZ_CC(Z_pos);
1432 pt2_recCluster->setRecZ_mTPC(tZ_pos);
1433 aCgemClusterCol->push_back(pt2_recCluster);
1434 // <--- fill 2D RecCgemCluster ---
1435
1436 if(nCluster[i][2]<100) {
1437 posZ[i][nCluster[i][2]]=Z_pos;
1438 tposZ[i][nCluster[i][2]]=tZ_pos;
1439 phi_XV[i][nCluster[i][2]]=vecPos_cluster[i][j][0][iX];
1440 }
1441 if( (iEnd_cluster[i][j][0][iX]+1)==readoutPlane->getNXstrips() )
1442 idxBoundaryXVcluster[i][j][0].push_back(idxCluster[i][j][1].size()-1);
1443 if(iStart_cluster[i][j][0][iX]==0)
1444 idxBoundaryXVcluster[i][j][1].push_back(idxCluster[i][j][1].size()-1);
1445 nCluster[i][2]++;
1446 }
1447 }// X-clusters loop ends
1448 // <--- loop X-V combinations ---
1449 // ------------------------------
1450 }
1451 nCluster[i][k]++;
1452
1453 // --- reset ---
1454 //resetForCluster();
1455 clusterFound=true;
1456 sumQ=0.0;
1457 sumQX=0.0;
1458 pos_strips.clear();
1459 drift_strips.clear();
1460 q_strips.clear();
1461 id_strips.clear();
1462 sum_x_tpc=0.;
1463 sum_y_tpc=0.;
1464 sum_xx_tpc=0.;
1465 sum_yy_tpc=0.;
1466 sum_xy_tpc=0.;
1467 n_tpc=0;
1468 a_tpc=0.;
1469 b_tpc=0.;
1470 pos_tpc=9999.;
1471 // -------------
1472 }
1473
1474 }// loop fired strips
1475 }// loop XV
1476 }// loop sheets
1477
1478 // to combine clusters near the boundary between sheets
1479 for(int j=0; j<myNSheets[i]; j++)
1480 {
1481 // get the index of the next sheet
1482 int j_next=j+1;
1483 if(j_next==myNSheets[i]) j_next=0;
1484
1485 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(i,j);
1486 CgemGeoReadoutPlane* nextReadoutPlane = myCgemGeomSvc->getReadoutPlane(i,j_next);
1487 double xmin_next = nextReadoutPlane->getXmin();
1488
1489 int nXV_atXEnd=idxBoundaryXVcluster[i][j][0].size();
1490 int nXV_atXStart=idxBoundaryXVcluster[i][j_next][1].size();
1491 if(nXV_atXEnd==0||nXV_atXStart==0) continue;
1492 for(int iXV_atXEnd=0; iXV_atXEnd<nXV_atXEnd; iXV_atXEnd++)
1493 {
1494 int iXVCluster = idxBoundaryXVcluster[i][j][0][iXV_atXEnd];
1495 int iXCluster = idxCluster[i][j][0][iXVCluster];
1496 int iVCluster = idxCluster[i][j][1][iXVCluster];
1497 int VID1=iStart_cluster[i][j][1][iVCluster];
1498 int VID2=iEnd_cluster[i][j][1][iVCluster];
1499 int VID1_extend=readoutPlane->getVIDInNextSheetFromVID(VID1, xmin_next);
1500 int VID2_extend=readoutPlane->getVIDInNextSheetFromVID(VID2, xmin_next);
1501 for(int iXV_atXStart=0; iXV_atXStart<nXV_atXStart; iXV_atXStart++)
1502 {
1503 int iXVCluster_next = idxBoundaryXVcluster[i][j_next][1][iXV_atXStart];
1504 int iXCluster_next = idxCluster[i][j_next][0][iXVCluster_next];
1505 int iVCluster_next = idxCluster[i][j_next][1][iXVCluster_next];
1506 int VID1_next=iStart_cluster[i][j_next][1][iVCluster_next];
1507 int VID2_next=iEnd_cluster[i][j_next][1][iVCluster_next];
1508 if( !( (VID1_next-VID2_extend)>1 || (VID1_extend-VID2_next)>1 ) )
1509 {
1510 int XID1=iStart_cluster[i][j][0][iXCluster];
1511 int XID2=iEnd_cluster[i][j][0][iXCluster];
1512 int XID1_next=iStart_cluster[i][j_next][0][iXCluster_next];
1513 int XID2_next=iEnd_cluster[i][j_next][0][iXCluster_next];
1514 int clusterId=aCgemClusterCol->size();
1515
1516 RecCgemCluster *pt_recCluster = new RecCgemCluster();
1517 pt_recCluster->setclusterid(clusterId);
1518 pt_recCluster->setlayerid(i);
1519 pt_recCluster->setsheetid(-1);
1520 pt_recCluster->setflag(3);
1521 int id1=id_cluster[i][j][2][iXVCluster];
1522 int id2=id_cluster[i][j_next][2][iXVCluster_next];
1523 pt_recCluster->setclusterflag(id1,id2);
1524 // averaging phi
1525 double phi1=vecPos_cluster[i][j][0][iXCluster];
1526 double phi2=vecPos_cluster[i][j_next][0][iXCluster_next];
1527 if(fabs(phi1-phi2)>CLHEP::pi) // phi domain 0~2pi, so only for 0/2pi
1528 {
1529 if(phi1>CLHEP::pi) phi1=phi1-CLHEP::twopi;
1530 if(phi2>CLHEP::pi) phi2=phi2-CLHEP::twopi;
1531 }
1532 double E1=E_cluster[i][j][0][iXCluster];
1533 double E2=E_cluster[i][j_next][0][iXCluster_next];
1534 double phiAverage=(E1*phi1+E2*phi2)/(E1+E2);
1535 pt_recCluster->setrecphi(phiAverage);
1536 // averging V
1537 double v1=vecPos_cluster[i][j][1][iVCluster];
1538 v1=readoutPlane->getVInNextSheetFromV(v1, xmin_next);
1539 double v2=vecPos_cluster[i][j_next][1][iVCluster_next];
1540 E1=E_cluster[i][j][1][iVCluster];
1541 E2=E_cluster[i][j_next][1][iVCluster_next];
1542 double V_average_next=(E1*v1+E2*v2)/(E1+E2);
1543 pt_recCluster->setrecv(V_average_next);
1544 // get Z
1545 double z_average = nextReadoutPlane->getZFromPhiV(phiAverage,V_average_next,0);
1546 pt_recCluster->setRecZ(z_average);
1547 // push back into CgemClusterCol
1548 aCgemClusterCol->push_back(pt_recCluster);
1549
1550 if(myPrintFlag) {
1551 cout<<"one combinational cluster found: "<<endl;
1552 cout<<" sheet "<<j <<" xID:"<<XID1 <<"~"<<XID2 <<", phi:"<<vecPos_cluster[i][j][0][iXCluster] <<", vID:"<<VID1 <<"~"<<VID2 <<", v:"<<vecPos_cluster[i][j][1][iVCluster] <<", z:"<<vecPos_cluster[i][j][2][iXVCluster] <<endl;
1553 cout<<" sheet "<<j_next<<" xID:"<<XID1_next<<"~"<<XID2_next<<", phi:"<<vecPos_cluster[i][j_next][0][iXCluster_next]<<", vID:"<<VID1_next<<"~"<<VID2_next<<", v:"<<vecPos_cluster[i][j_next][1][iVCluster_next]<<", z:"<<vecPos_cluster[i][j_next][2][iXVCluster_next]<<endl;
1554 cout<<" averagePhi:"<<phiAverage<<", v_average:"<<V_average_next<<", z_average:"<<z_average<<endl;
1555 }
1556 }
1557 }
1558 }
1559 }// end combining clusters near the boundary between sheets
1560
1561 if(nCluster[i][0]>100) nCluster[i][0]=100;
1562 if(nCluster[i][1]>100) nCluster[i][1]=100;
1563 if(nCluster[i][2]>100) nCluster[i][2]=100;
1564 }// loop layers
1565
1566 // count time used
1567 m_timer->stop();
1568 double evttime=m_timer->elapsed();
1569
1570 // checkRecCgemClusterCol
1571 if(myPrintFlag) checkRecCgemClusterCol();
1572
1573 if(myNtuple)
1574 {
1575
1576 myNTupleHelper->fillLong("run",(long)run);
1577 myNTupleHelper->fillLong("evt",(long)evt);
1578
1579 myNTupleHelper->fillArray("nXstrips","nLayer",(double*) nXStrips,3);
1580 myNTupleHelper->fillArray("nVstrips","nLayer",(double*) nVStrips,3);
1581
1582 myNTupleHelper->fillArray("phiClusterLay1","nXClusterLay1",(double*) posX[0],nCluster[0][0]);
1583 myNTupleHelper->fillArray("tphiClusterLay1","nXClusterLay1",(double*) tposX[0],nCluster[0][0]);
1584 myNTupleHelper->fillArrayInt("nXStripsLay1","nXClusterLay1",(int*) nStripsX[0],nCluster[0][0]);
1585 myNTupleHelper->fillArray("QXLay1","nXClusterLay1",(double*) QX[0],nCluster[0][0]);
1586 myNTupleHelper->fillArray("a_tpc_XLay1" ,"nXClusterLay1",(double*) a_tpc_cluster_X[0] ,nCluster[0][0]);
1587 myNTupleHelper->fillArray("b_tpc_XLay1" ,"nXClusterLay1",(double*) b_tpc_cluster_X[0] ,nCluster[0][0]);
1588 myNTupleHelper->fillArray("pos_tpc_XLay1" ,"nXClusterLay1",(double*) pos_tpc_cluster_X[0],nCluster[0][0]);
1589
1590 myNTupleHelper->fillArray("VClusterLay1","nVClusterLay1",(double*) posV[0],nCluster[0][1]);
1591 myNTupleHelper->fillArray("tVClusterLay1","nVClusterLay1",(double*) tposV[0],nCluster[0][1]);
1592 myNTupleHelper->fillArrayInt("nVStripsLay1","nVClusterLay1",(int*) nStripsV[0],nCluster[0][1]);
1593 myNTupleHelper->fillArray("QVLay1","nVClusterLay1",(double*) QV[0],nCluster[0][1]);
1594
1595 myNTupleHelper->fillArray("zClusterLay1","nXVClusterLay1",(double*) posZ[0],nCluster[0][2]);
1596 myNTupleHelper->fillArray("tzClusterLay1","nXVClusterLay1",(double*) tposZ[0],nCluster[0][2]);
1597 myNTupleHelper->fillArray("phiXVClusterLay1","nXVClusterLay1",(double*) phi_XV[0],nCluster[0][2]);
1598
1599 myNTupleHelper->fillArray("phiClusterLay2","nXClusterLay2",(double*) posX[1],nCluster[1][0]);
1600 myNTupleHelper->fillArray("tphiClusterLay2","nXClusterLay2",(double*) tposX[1],nCluster[1][0]);
1601 myNTupleHelper->fillArrayInt("nXStripsLay2","nXClusterLay2",(int*) nStripsX[1],nCluster[1][0]);
1602 myNTupleHelper->fillArray("QXLay2","nXClusterLay2",(double*) QX[1],nCluster[1][0]);
1603 myNTupleHelper->fillArray("a_tpc_XLay2" ,"nXClusterLay2",(double*) a_tpc_cluster_X[1] ,nCluster[1][0]);
1604 myNTupleHelper->fillArray("b_tpc_XLay2" ,"nXClusterLay2",(double*) b_tpc_cluster_X[1] ,nCluster[1][0]);
1605 myNTupleHelper->fillArray("pos_tpc_XLay2" ,"nXClusterLay2",(double*) pos_tpc_cluster_X[1],nCluster[1][0]);
1606
1607 myNTupleHelper->fillArray("VClusterLay2","nVClusterLay2",(double*) posV[1],nCluster[1][1]);
1608 myNTupleHelper->fillArray("tVClusterLay2","nVClusterLay2",(double*) tposV[1],nCluster[1][1]);
1609 myNTupleHelper->fillArrayInt("nVStripsLay2","nVClusterLay2",(int*) nStripsV[1],nCluster[1][1]);
1610 myNTupleHelper->fillArray("QVLay2","nVClusterLay2",(double*) QV[1],nCluster[1][1]);
1611
1612 myNTupleHelper->fillArray("zClusterLay2","nXVClusterLay2",(double*) posZ[1],nCluster[1][2]);
1613 myNTupleHelper->fillArray("tzClusterLay2","nXVClusterLay2",(double*) tposZ[1],nCluster[1][2]);
1614 myNTupleHelper->fillArray("phiXVClusterLay2","nXVClusterLay2",(double*) phi_XV[1],nCluster[1][2]);
1615
1616 myNTupleHelper->fillArray("phiClusterLay3","nXClusterLay3",(double*) posX[2],nCluster[2][0]);
1617 myNTupleHelper->fillArray("tphiClusterLay3","nXClusterLay3",(double*) tposX[2],nCluster[2][0]);
1618 myNTupleHelper->fillArrayInt("nXStripsLay3","nXClusterLay3",(int*) nStripsX[2],nCluster[2][0]);
1619 myNTupleHelper->fillArray("QXLay3","nXClusterLay3",(double*) QX[2],nCluster[2][0]);
1620 myNTupleHelper->fillArray("a_tpc_XLay3" ,"nXClusterLay3",(double*) a_tpc_cluster_X[2] ,nCluster[2][0]);
1621 myNTupleHelper->fillArray("b_tpc_XLay3" ,"nXClusterLay3",(double*) b_tpc_cluster_X[2] ,nCluster[2][0]);
1622 myNTupleHelper->fillArray("pos_tpc_XLay3" ,"nXClusterLay3",(double*) pos_tpc_cluster_X[2],nCluster[2][0]);
1623
1624 myNTupleHelper->fillArray("VClusterLay3","nVClusterLay3",(double*) posV[2],nCluster[2][1]);
1625 myNTupleHelper->fillArray("tVClusterLay3","nVClusterLay3",(double*) tposV[2],nCluster[2][1]);
1626 myNTupleHelper->fillArrayInt("nVStripsLay3","nVClusterLay3",(int*) nStripsV[2],nCluster[2][1]);
1627 myNTupleHelper->fillArray("QVLay3","nVClusterLay3",(double*) QV[2],nCluster[2][1]);
1628
1629 myNTupleHelper->fillArray("zClusterLay3","nXVClusterLay3",(double*) posZ[2],nCluster[2][2]);
1630 myNTupleHelper->fillArray("tzClusterLay3","nXVClusterLay3",(double*) tposZ[2],nCluster[2][2]);
1631 myNTupleHelper->fillArray("phiXVClusterLay3","nXVClusterLay3",(double*) phi_XV[2],nCluster[2][2]);
1632
1633 myNTupleHelper->fillDouble("evtProcessTime",evttime);
1634
1635 myNTupleHelper->write();
1636 }
1637 if(myPrintFlag) cout<<"End of CgemClusterCreate::method2()"<<endl;
1638
1639 return StatusCode::SUCCESS;
1640}// method2
1641
1642StatusCode CgemClusterCreate::toyCluster()
1643{
1644 MsgStream log(msgSvc(),name());
1645
1646 // --- get evtHeader
1647 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
1648 if (!evtHead)
1649 {
1650 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
1651 return StatusCode::FAILURE;
1652 }
1653 int run=evtHead->runNumber();
1654 int evt=evtHead->eventNumber();
1655 if(myPrintFlag) cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ "<<evt<<endl;
1656
1657 // --- register a RecCgemClusterCol
1658 IDataProviderSvc* evtSvc = NULL;
1659 Gaudi::svcLocator()->service("EventDataSvc",evtSvc);
1660 if (evtSvc) {
1661 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1662 } else {
1663 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1664 return StatusCode::SUCCESS;
1665 }
1666 StatusCode sc;
1667 RecCgemClusterCol* aCgemClusterCol = new RecCgemClusterCol;
1668 sc = evtSvc->registerObject("/Event/Recon/RecCgemClusterCol", aCgemClusterCol);
1669 if(sc.isFailure())
1670 {
1671 log << MSG::FATAL << "Could not register RecCgemClusterCol" << endreq;
1672 return StatusCode::SUCCESS;
1673 }
1674
1675 // --- get CgemMcHitCol
1676 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(),"/Event/MC/CgemMcHitCol");
1677 if (!cgemMcHitCol)
1678 {
1679 log << MSG::WARNING << "Could not retrieve CgemMcHitCol list" << endreq;
1680 return StatusCode::FAILURE;
1681 }
1682
1683
1684 // --- loop CgemMcHit and generate CgemCluster
1685 double phi_truth[100], z_truth[100], r_truth[100], x_truth[100], v_truth[100], cosTh_truth[100], z_check[100];
1686 double phi_gen[100], z_gen[100], x_gen[100], v_gen[100];
1687 int iCgemMcHit(0);
1688 int nCgemMcHit = cgemMcHitCol->size(); if(myPrintFlag) cout<<"nCgemMcHit = "<<nCgemMcHit<<endl;
1689 vector<int> idxXClusters[4][2], idxVClusters[4][2];
1690 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
1691
1692 for(; iter_truth!=cgemMcHitCol->end(); iter_truth++ )
1693 {
1694 if(myPrintFlag) cout<<"CgemClusterCreate::toyCluster(): creator process = "<<(*iter_truth)->GetCreatorProcess()<<endl;
1695 string creatorProcess = (*iter_truth)->GetCreatorProcess();
1696 if(creatorProcess=="Generator"||creatorProcess=="Decay")
1697 {
1698 // sampling the efficiency of clustering
1699 if(myClusterEff>0.&&myClusterEff<1.0) {
1700 Rndm::Numbers flat(randSvc(), Rndm::Flat(0,1));
1701 if(flat()>myClusterEff) {
1702 if(myPrintFlag) cout<<"CgemClusterCreate::toyCluster() skip one cluster! "<<endl;
1703 continue;
1704 }
1705 }
1706 /* get position of CgemMcHit (in mm)*/
1707 int iLayer = (*iter_truth)->GetLayerID();
1708 double x_pre = (*iter_truth)->GetPositionXOfPrePoint();
1709 double y_pre = (*iter_truth)->GetPositionYOfPrePoint();
1710 double z_pre = (*iter_truth)->GetPositionZOfPrePoint();
1711 double x_post = (*iter_truth)->GetPositionXOfPostPoint();
1712 double y_post = (*iter_truth)->GetPositionYOfPostPoint();
1713 double z_post = (*iter_truth)->GetPositionZOfPostPoint();
1714 double x_mid = 0.5*(x_pre+x_post);
1715 double y_mid = 0.5*(y_pre+y_post);
1716 double z_mid = 0.5*(z_pre+z_post);
1717 double r_pre = sqrt(x_pre*x_pre+y_pre*y_pre+z_pre*z_pre);
1718 double r_post = sqrt(x_post*x_post+y_post*y_post+z_post*z_post);
1719 double v_x = x_post-x_pre;
1720 double v_y = y_post-y_pre;
1721 double v_z = z_post-z_pre;
1722 double cth_v = v_z/sqrt(v_x*v_x+v_y*v_y+v_z*v_z);
1723 if(r_pre>r_post) {
1724 v_x*=-1.0;
1725 v_y*=-1.0;
1726 v_z*=-1.0;
1727 }
1728 double v_tot = sqrt(v_x*v_x+v_y*v_y+v_z*v_z);
1729 double theta_v = acos(v_z/v_tot);
1730 double phi_v = atan2(v_y, v_x);
1731 //double r_middle = sqrt(0.25*(x_pre+x_post)*(x_pre+x_post)+0.25*(y_pre+y_post)*(y_pre+y_post));
1732 double r_middle = sqrt(x_mid*x_mid+y_mid*y_mid);
1733 if(myPrintFlag) cout<<"iLayer, r_middle = "<<iLayer<<", "<<r_middle<<endl;
1734 Hep3Vector pos_mid(x_mid, y_mid, z_mid);
1735 double phi_mid = pos_mid.phi();
1736 while(phi_mid>CLHEP::twopi) phi_mid-=CLHEP::twopi;
1737 while(phi_mid<0.) phi_mid+=CLHEP::twopi;
1738 /* angle between track and radial direction */
1739 double dPhi = phi_v-phi_mid;
1740 while(dPhi>CLHEP::pi) dPhi-=CLHEP::twopi;
1741 while(dPhi<-CLHEP::pi) dPhi+=CLHEP::twopi;
1742 /* get readout plane */
1743 CgemGeoReadoutPlane* readoutPlane=NULL;
1744 int iSheet=0;
1745 for(int i=0; i<myNSheets[iLayer]; i++)
1746 {
1747 readoutPlane = myCgemGeomSvc->getReadoutPlane(iLayer,i);
1748 if(readoutPlane->OnThePlane(phi_mid,z_mid)) // rad, mm
1749 {
1750 iSheet = i;
1751 break;
1752 }
1753 readoutPlane=NULL;
1754 }
1755 if(readoutPlane==NULL)continue;
1756 /* get V truth */
1757 double v_loc = readoutPlane->getVFromPhiZ(phi_mid, z_mid);// mm
1758 if(iCgemMcHit<100) {
1759 phi_truth[iCgemMcHit]=phi_mid;
1760 z_truth[iCgemMcHit]=z_mid;//in mm
1761 r_truth[iCgemMcHit]=r_middle;//in mm
1762 x_truth[iCgemMcHit]=phi_mid*myRXstrips[iLayer];//mm
1763 v_truth[iCgemMcHit]=v_loc;//mm
1764 cosTh_truth[iCgemMcHit]=cth_v;
1765 z_check[iCgemMcHit] = readoutPlane->getZFromPhiV(phi_mid, v_loc);
1766 }
1767 if(myPrintFlag) cout<<"CgemClusterCreate::toyCluster() iLayer "<<iLayer<<", MC hit: phi, z, V = "<<phi_mid<<", "<<z_mid<<", "<<v_loc<<endl;
1768
1769 /* generate new position */
1770 // phi view
1771 int iView(0), mode(2);
1772 double Q(100), T(100);
1773 double sigma_X = myCgemCalibSvc->getSigma(iLayer,iView,mode,dPhi,Q,T);// in mm
1774 Rndm::Numbers gauss(randSvc(), Rndm::Gauss(0,1));
1775 double phi_mid_gen = phi_mid+gauss()*sigma_X/myRXstrips[iLayer];// in radian
1776 int iGenTried=0;
1777 while(!(readoutPlane->OnThePlane(phi_mid_gen,z_mid))) {
1778 phi_mid_gen = phi_mid+gauss()*sigma_X/myRXstrips[iLayer];
1779 iGenTried++;
1780 if(iGenTried>=10) break;
1781 }
1782 if(iGenTried>=10) {
1783 if(myPrintFlag) cout<<"Generation of phi with 10 times!"<<endl;
1784 continue;
1785 }
1786 if(myPrintFlag) cout<<"generated phi: "<<phi_mid_gen<<endl;
1787 // V view
1788 iView=1;
1789 double sigma_V = myCgemCalibSvc->getSigma(iLayer,iView,mode,dPhi,Q,T);// in mm
1790 //sigma_V = sigma_V/10.;// in cm
1791 double v_loc_gen = v_loc+gauss()*sigma_V;// mm
1792 double z_mid_gen = readoutPlane->getZFromPhiV(phi_mid_gen, v_loc_gen);
1793 iGenTried=0;
1794 while(!(readoutPlane->OnThePlane(phi_mid_gen,z_mid_gen)))
1795 {
1796 v_loc_gen = v_loc+gauss()*sigma_V;
1797 z_mid_gen = readoutPlane->getZFromPhiV(phi_mid_gen, v_loc_gen);
1798 if(myPrintFlag) cout<<"try generated V, z: "<<v_loc_gen<<", "<<z_mid_gen<<endl;
1799 iGenTried++;
1800 if(iGenTried>=10) break;
1801 }
1802 if(iGenTried>=10) {
1803 if(myPrintFlag) cout<<"Generation of V with 10 times!"<<endl;
1804 continue;
1805 }
1806 //z_mid = z_mid+gauss()*sigma_V;// in mm
1807 if(iCgemMcHit<100) {
1808 phi_gen[iCgemMcHit]=phi_mid_gen;
1809 z_gen[iCgemMcHit]=z_mid_gen;
1810 x_gen[iCgemMcHit]=phi_mid_gen*myRXstrips[iLayer];
1811 v_gen[iCgemMcHit]=v_loc_gen;
1812
1813 }
1814
1815 /** fill RecCgemCluster */
1816 /* fill X */
1817 int clusterId=aCgemClusterCol->size();
1818 RecCgemCluster *pt_recCluster = new RecCgemCluster();
1819 //cout<<"clusterId = "<<clusterId<<endl;
1820 pt_recCluster->setclusterid(clusterId);
1821 pt_recCluster->setlayerid(iLayer);
1822 pt_recCluster->setsheetid(iSheet);
1823 pt_recCluster->setflag(0);
1824 //pt_recCluster->setclusterflag(id_cluster[i][j][0][iX],id_cluster[i][j][1][iV]);
1825 pt_recCluster->setrecphi(phi_mid_gen);
1826 //pt_recCluster->setrecv(vecPos_cluster[i][j][1][iV]);
1827 //pt_recCluster->setRecZ(z_mid*10);// in mm
1828 aCgemClusterCol->push_back(pt_recCluster);
1829 idxXClusters[iLayer][iSheet].push_back(clusterId);
1830 (*iter_truth)->AddXclusterIdx(clusterId);// associate X-cluster id to CGEM true hit
1831 /* fill V */
1832 clusterId=aCgemClusterCol->size();
1833 pt_recCluster = new RecCgemCluster();
1834 //cout<<"clusterId = "<<clusterId<<endl;
1835 pt_recCluster->setclusterid(clusterId);
1836 pt_recCluster->setlayerid(iLayer);
1837 pt_recCluster->setsheetid(iSheet);
1838 pt_recCluster->setflag(1);
1839 pt_recCluster->setrecv(v_loc_gen);// mm
1840 aCgemClusterCol->push_back(pt_recCluster);
1841 idxVClusters[iLayer][iSheet].push_back(clusterId);
1842 (*iter_truth)->AddVclusterIdx(clusterId);// associate V-cluster id to CGEM true hit
1843
1844 // --- momentum
1845 int pdg = (*iter_truth)->GetPDGCode();
1846 double px = (*iter_truth)->GetMomentumXOfPrePoint();
1847 double py = (*iter_truth)->GetMomentumYOfPrePoint();
1848 double pz = (*iter_truth)->GetMomentumZOfPrePoint();
1849 Hep3Vector truth_p(px/1000.,py/1000.,pz/1000.);
1850 //cout<<"pt="<<sqrt(px*px+py*py)<<endl;
1851 HepPoint3D truth_position(x_pre/10.,y_pre/10.,z_pre/10.);
1852 if(myPrintFlag&&fabs(pdg)==13)
1853 {
1854 int charge = -1*pdg/fabs(pdg);
1855 KalmanFit::Helix* tmp_helix = new KalmanFit::Helix(truth_position, truth_p, charge);
1856 HepPoint3D tmp_pivot(0., 0., 0.);
1857 tmp_helix->pivot(tmp_pivot);
1858 if(myPrintFlag)
1859 cout<<"pdg="<<pdg<<setw(10)<<" Helix of MC truth: "<<setw(10)<<tmp_helix->dr()<<setw(10)<<tmp_helix->phi0()<<setw(10)<<tmp_helix->kappa()<<setw(10)<<tmp_helix->dz()<<setw(10)<<tmp_helix->tanl()<<endl;
1860 delete tmp_helix;
1861 }
1862
1863 iCgemMcHit++;
1864
1865 } // if creatorProcess
1866 } // loop CgemMcHit
1867
1868 // --- add random noise
1869 if(myRandomCluster>0) {
1870 Rndm::Numbers flat(randSvc(), Rndm::Flat(0,1));
1871 for(int i=0; i<myNCgemLayers; i++) {// loop layers
1872 Rndm::Numbers poisson(randSvc(), Rndm::Poisson(myNoiseRate[i]/myNSheets[i]));
1873 for(int j=0; j<myNSheets[i]; j++)// loop sheet to generate random 1D clusters
1874 {
1875 CgemGeoReadoutPlane* readoutPlane = myCgemGeomSvc->getReadoutPlane(i,j);
1876 int nBK=poisson();
1877 for(int k=0; k<nBK; k++) {
1878 double phi=flat()*(2.0*M_PI)/myNSheets[i]+(j-1)*M_PI;
1879 double z=flat()*readoutPlane->getLength()+readoutPlane->getZmin();
1880 if(readoutPlane->OnThePlane(phi,z)) // rad, mm
1881 {
1882 /** fill RecCgemCluster */
1883 /* fill X */
1884 int clusterId=aCgemClusterCol->size();
1885 RecCgemCluster *pt_recCluster = new RecCgemCluster();
1886 //cout<<"clusterId = "<<clusterId<<endl;
1887 pt_recCluster->setclusterid(clusterId);
1888 pt_recCluster->setlayerid(i);
1889 pt_recCluster->setsheetid(j);
1890 pt_recCluster->setflag(0);
1891 pt_recCluster->setrecphi(phi);
1892 aCgemClusterCol->push_back(pt_recCluster);
1893 idxXClusters[i][j].push_back(clusterId);
1894
1895 /* fill V */
1896 double v_loc = readoutPlane->getVFromPhiZ(phi, z);// mm
1897 clusterId=aCgemClusterCol->size();
1898 pt_recCluster = new RecCgemCluster();
1899 //cout<<"clusterId = "<<clusterId<<endl;
1900 pt_recCluster->setclusterid(clusterId);
1901 pt_recCluster->setlayerid(i);
1902 pt_recCluster->setsheetid(j);
1903 pt_recCluster->setflag(1);
1904 pt_recCluster->setrecv(v_loc);// mm
1905 aCgemClusterCol->push_back(pt_recCluster);
1906 idxVClusters[i][j].push_back(clusterId);
1907 }
1908 }
1909 }
1910 }// loop layers
1911 }
1912
1913 //checkRecCgemClusterCol();
1914 /* search XV 2D clusters */
1915 int nXVCluster(0);
1916 RecCgemClusterCol::iterator it_cluster0 = aCgemClusterCol->begin();
1917 RecCgemClusterCol::iterator it_cluster;
1918 if(myPrintFlag) cout<<"CgemClusterCreate::toyCluster() start searching XV clusters: "<<endl;
1919 for(int i=0; i<myNCgemLayers; i++)
1920 {
1921 for(int j=0; j<myNSheets[i]; j++)
1922 {
1923 if(myPrintFlag) cout<<"iLayer, iSheet = "<<i<<", "<<j<<endl;
1924 CgemGeoReadoutPlane* readoutPlane= myCgemGeomSvc->getReadoutPlane(i,j);
1925 for(int iV=0; iV<idxVClusters[i][j].size(); iV++)
1926 {
1927 if(myPrintFlag) cout<<"iV: "<<iV;
1928 it_cluster = it_cluster0+idxVClusters[i][j][iV];
1929 //cout<<"get it_cluster"<<endl;
1930 double V_loc = (*it_cluster)->getrecv();
1931 //cout<<"get V_loc"<<endl;
1932 for(int iX=0; iX<idxXClusters[i][j].size(); iX++)
1933 {
1934 if(myPrintFlag) cout<<" iX: "<<iX<<endl;
1935 it_cluster = it_cluster0+idxXClusters[i][j][iX];
1936 double phi = (*it_cluster)->getrecphi();
1937 double z = readoutPlane->getZFromPhiV(phi,V_loc);
1938 if(myPrintFlag) cout<<"CgemClusterCreate::toyCluster() check phi, z, V = "<<phi<<", "<<z<<", "<<V_loc<<", layer "<<i<<", sheet "<<j<<endl;
1939 if(readoutPlane->OnThePlane(phi,z)) // 2018-05-15 by llwang
1940 {
1941 int clusterId=aCgemClusterCol->size();
1942 RecCgemCluster *pt_recCluster = new RecCgemCluster();
1943 //cout<<"clusterId = "<<clusterId<<endl;
1944 pt_recCluster->setclusterid(clusterId);
1945 pt_recCluster->setlayerid(i);
1946 pt_recCluster->setsheetid(j);
1947 pt_recCluster->setflag(2);
1948 pt_recCluster->setclusterflag(idxXClusters[i][j][iX], idxVClusters[i][j][iV]);
1949 pt_recCluster->setrecphi(phi);
1950 pt_recCluster->setrecv(V_loc);
1951 pt_recCluster->setRecZ(z);// in mm
1952 aCgemClusterCol->push_back(pt_recCluster);
1953 it_cluster0 = aCgemClusterCol->begin();// update it_cluster0
1954 if(myPrintFlag) cout<<"CgemClusterCreate::toyCluster() finds a XV cluster"<<endl;
1955 nXVCluster++;
1956 }
1957 }
1958 }
1959 }
1960 }// make 2D clusters
1961
1962 if(myPrintFlag) cout<<"nCgemCluster = "<<aCgemClusterCol->size()<<endl;
1963 //RecCgemClusterCol::iterator it_Cluster = aCgemClusterCol-
1964 if(myPrintFlag) checkRecCgemClusterCol();
1965 if(run<0) fillMCTruth();
1966 if(myNtuple)
1967 {
1968 if(iCgemMcHit>100) iCgemMcHit=100;
1969 myNTupleHelper->fillArray("phi_mc","nCgemMcHit",(double*)phi_truth,iCgemMcHit);
1970 myNTupleHelper->fillArray("z_mc","nCgemMcHit",(double*)z_truth,iCgemMcHit);
1971 myNTupleHelper->fillArray("z_mc_check","nCgemMcHit",(double*)z_check,iCgemMcHit);
1972 myNTupleHelper->fillArray("r_mc","nCgemMcHit",(double*)r_truth,iCgemMcHit);
1973 myNTupleHelper->fillArray("x_mc","nCgemMcHit",(double*)x_truth,iCgemMcHit);
1974 myNTupleHelper->fillArray("v_mc","nCgemMcHit",(double*)v_truth,iCgemMcHit);
1975 myNTupleHelper->fillArray("cth_mc","nCgemMcHit",(double*)cosTh_truth,iCgemMcHit);
1976 myNTupleHelper->fillArray("phi","nCgemMcHit",(double*)phi_gen,iCgemMcHit);
1977 myNTupleHelper->fillArray("z","nCgemMcHit",(double*)z_gen,iCgemMcHit);
1978 myNTupleHelper->fillArray("x","nCgemMcHit",(double*)x_gen,iCgemMcHit);
1979 myNTupleHelper->fillArray("v","nCgemMcHit",(double*)v_gen,iCgemMcHit);
1980 myNTupleHelper->write();
1981 //cout<<"CgemClusterCreate::myNTupleHelper filled!"<<endl;
1982 }//
1983 /*
1984 if(nXVCluster ==3){
1985 if(myNtuple)fillRecData(run,evt);
1986 if(myNtuple&&run<0)fillMCTruth(run,evt);
1987 }else
1988 cout<<"event:"<<evt<<" nXVCluster = "<<nXVCluster<<endl;
1989 */
1990 return StatusCode::SUCCESS;
1991}
1992
1993void CgemClusterCreate::fillMCTruth(int run, int evt)
1994{
1995 m_mc_run = run;
1996 m_mc_evt = evt;
1997 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
1998 if (!cgemMcHitCol)
1999 {
2000 std::cout << "Could not retrieve cgemMcHitCol" << std::endl;
2001 return;
2002 }else{
2003 m_mc_nhit = cgemMcHitCol->size();
2004 //if(m_mc_nhit>3)cout<<evt<<" "<<m_mc_nhit<<endl;
2005 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
2006 for(int iHit=0; iter_truth!= cgemMcHitCol->end(); iter_truth++,iHit++)
2007 {
2008 m_mc_trackID[iHit] = (*iter_truth)->GetTrackID();
2009 m_mc_layerID[iHit] = (*iter_truth)->GetLayerID();
2010 m_mc_pdg[iHit] = (*iter_truth)->GetPDGCode();
2011 m_mc_parentID[iHit] = (*iter_truth)->GetParentID();
2012 m_mc_E_deposit[iHit] = (*iter_truth)->GetTotalEnergyDeposit();
2013 m_mc_XYZ_pre_x[iHit] = (*iter_truth)->GetPositionXOfPrePoint();
2014 m_mc_XYZ_pre_y[iHit] = (*iter_truth)->GetPositionYOfPrePoint();
2015 m_mc_XYZ_pre_z[iHit] = (*iter_truth)->GetPositionZOfPrePoint();
2016 m_mc_XYZ_post_x[iHit] = (*iter_truth)->GetPositionXOfPostPoint();
2017 m_mc_XYZ_post_y[iHit] = (*iter_truth)->GetPositionYOfPostPoint();
2018 m_mc_XYZ_post_z[iHit] = (*iter_truth)->GetPositionZOfPostPoint();
2019 m_mc_P_pre_x[iHit] = (*iter_truth)->GetMomentumXOfPrePoint();
2020 m_mc_P_pre_y[iHit] = (*iter_truth)->GetMomentumYOfPrePoint();
2021 m_mc_P_pre_z[iHit] = (*iter_truth)->GetMomentumZOfPrePoint();
2022 m_mc_P_post_x[iHit] = (*iter_truth)->GetMomentumXOfPostPoint();
2023 m_mc_P_post_y[iHit] = (*iter_truth)->GetMomentumYOfPostPoint();
2024 m_mc_P_post_z[iHit] = (*iter_truth)->GetMomentumZOfPostPoint();
2025 if(iHit>1000)break;
2026 }
2027 }
2028 m_mc_nt->write();
2029}
2030
2031void CgemClusterCreate::fillRecData(int run, int evt)
2032{
2033 int iCluster(0);
2034 //int fillOption = -1;
2035 m_rec_run = run;
2036 m_rec_evt = evt;
2037 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
2038 if(!recCgemClusterCol)
2039 {
2040 cout << "Could not retrieve RecCgemClusterCol" << endl;
2041 }else{
2042 m_rec_ncluster = recCgemClusterCol->size();
2043 //if(m_rec_ncluster>3)cout<<m_rec_ncluster<<" ";
2044 RecCgemClusterCol::iterator iter_cluster=recCgemClusterCol->begin();
2045 for(; iter_cluster!=recCgemClusterCol->end(); iter_cluster++)
2046 {
2047 if(m_fillOption == -1){
2048 m_rec_clusterid[iCluster] = (*iter_cluster)->getclusterid();
2049 m_rec_layerid[iCluster] = (*iter_cluster)->getlayerid();
2050 m_rec_sheetid[iCluster] = (*iter_cluster)->getsheetid();
2051 m_rec_flag[iCluster] = (*iter_cluster)->getflag();
2052 m_rec_energydeposit[iCluster] = (*iter_cluster)->getenergydeposit();
2053 m_rec_recphi[iCluster] = (*iter_cluster)->getrecphi();
2054 m_rec_recv[iCluster] = (*iter_cluster)->getrecv();
2055 m_rec_recZ[iCluster] = (*iter_cluster)->getRecZ();
2056 m_rec_clusterflagb[iCluster] = (*iter_cluster)->getclusterflagb();
2057 m_rec_clusterflage[iCluster] = (*iter_cluster)->getclusterflage();
2058 iCluster++;
2059 }else if(m_fillOption==(*iter_cluster)->getflag()){
2060 m_rec_clusterid[iCluster] = (*iter_cluster)->getclusterid();
2061 m_rec_layerid[iCluster] = (*iter_cluster)->getlayerid();
2062 m_rec_sheetid[iCluster] = (*iter_cluster)->getsheetid();
2063 m_rec_flag[iCluster] = (*iter_cluster)->getflag();
2064 m_rec_energydeposit[iCluster] = (*iter_cluster)->getenergydeposit();
2065 m_rec_recphi[iCluster] = (*iter_cluster)->getrecphi();
2066 m_rec_recv[iCluster] = (*iter_cluster)->getrecv();
2067 m_rec_recZ[iCluster] = (*iter_cluster)->getRecZ();
2068 m_rec_clusterflagb[iCluster] = (*iter_cluster)->getclusterflagb();
2069 m_rec_clusterflage[iCluster] = (*iter_cluster)->getclusterflage();
2070 iCluster++;
2071 }
2072 if(iCluster>1000)break;
2073 }
2074 }
2075 m_rec_ncluster = iCluster;
2076 //if(iCluster>3)cout<<iCluster<<" ";
2077 m_rec_nt->write();
2078}
2079
2080void CgemClusterCreate::hist_def(){
2081 StatusCode status;
2082 NTuplePtr nt_rec(ntupleSvc(),"RecCgem/RecCluster");
2083 if( nt_rec ) m_rec_nt = nt_rec;
2084 else{
2085 m_rec_nt= ntupleSvc()->book("RecCgem/RecCluster",CLID_ColumnWiseTuple,"RecCluster");
2086 if( m_rec_nt ){
2087 status = m_rec_nt->addItem("rec_run" ,m_rec_run);
2088 status = m_rec_nt->addItem("rec_evt" ,m_rec_evt);
2089 status = m_rec_nt->addItem("rec_ncluster" ,m_rec_ncluster,0,1000);
2090 status = m_rec_nt->addItem("rec_clusterid" ,m_rec_ncluster,m_rec_clusterid);
2091 status = m_rec_nt->addItem("rec_layerid" ,m_rec_ncluster,m_rec_layerid);
2092 status = m_rec_nt->addItem("rec_sheetid" ,m_rec_ncluster,m_rec_sheetid);
2093 status = m_rec_nt->addItem("rec_flag" ,m_rec_ncluster,m_rec_flag);
2094 status = m_rec_nt->addItem("rec_energydeposit" ,m_rec_ncluster,m_rec_energydeposit);
2095 status = m_rec_nt->addItem("rec_recphi" ,m_rec_ncluster,m_rec_recphi);
2096 status = m_rec_nt->addItem("rec_recv" ,m_rec_ncluster,m_rec_recv);
2097 status = m_rec_nt->addItem("rec_recZ" ,m_rec_ncluster,m_rec_recZ);
2098 status = m_rec_nt->addItem("rec_clusterflagb" ,m_rec_ncluster,m_rec_clusterflagb);
2099 status = m_rec_nt->addItem("rec_clusterflage" ,m_rec_ncluster,m_rec_clusterflage);
2100 }
2101 }
2102
2103 NTuplePtr nt_mc(ntupleSvc(),"RecCgem/McHit");
2104 if( nt_mc ) m_mc_nt = nt_mc;
2105 else{
2106 m_mc_nt = ntupleSvc()->book("RecCgem/McHit",CLID_ColumnWiseTuple,"McHit");
2107 if( m_mc_nt ){
2108 status = m_mc_nt->addItem("mc_run" ,m_mc_run);
2109 status = m_mc_nt->addItem("mc_evt" ,m_mc_evt);
2110 status = m_mc_nt->addItem("mc_nhit" ,m_mc_nhit,0,1000);
2111 status = m_mc_nt->addItem("mc_trackID" ,m_mc_nhit,m_mc_trackID);
2112 status = m_mc_nt->addItem("mc_layerID" ,m_mc_nhit,m_mc_layerID);
2113 status = m_mc_nt->addItem("mc_pdg" ,m_mc_nhit,m_mc_pdg);
2114 status = m_mc_nt->addItem("mc_parentID" ,m_mc_nhit,m_mc_parentID);
2115 status = m_mc_nt->addItem("mc_E_deposit" ,m_mc_nhit,m_mc_E_deposit);
2116 status = m_mc_nt->addItem("mc_XYZ_pre_x" ,m_mc_nhit,m_mc_XYZ_pre_x);
2117 status = m_mc_nt->addItem("mc_XYZ_pre_y" ,m_mc_nhit,m_mc_XYZ_pre_y);
2118 status = m_mc_nt->addItem("mc_XYZ_pre_z" ,m_mc_nhit,m_mc_XYZ_pre_z);
2119 status = m_mc_nt->addItem("mc_XYZ_post_x" ,m_mc_nhit,m_mc_XYZ_post_x);
2120 status = m_mc_nt->addItem("mc_XYZ_post_y" ,m_mc_nhit,m_mc_XYZ_post_y);
2121 status = m_mc_nt->addItem("mc_XYZ_post_z" ,m_mc_nhit,m_mc_XYZ_post_z);
2122 status = m_mc_nt->addItem("mc_P_pre_x" ,m_mc_nhit,m_mc_P_pre_x);
2123 status = m_mc_nt->addItem("mc_P_pre_y" ,m_mc_nhit,m_mc_P_pre_y);
2124 status = m_mc_nt->addItem("mc_P_pre_z" ,m_mc_nhit,m_mc_P_pre_z);
2125 status = m_mc_nt->addItem("mc_P_post_x" ,m_mc_nhit,m_mc_P_post_x);
2126 status = m_mc_nt->addItem("mc_P_post_y" ,m_mc_nhit,m_mc_P_post_y);
2127 status = m_mc_nt->addItem("mc_P_post_z" ,m_mc_nhit,m_mc_P_post_z);
2128 }
2129 }
2130}
2131
2132StatusCode CgemClusterCreate::finalize()
2133{
2134 MsgStream log(msgSvc(),name());
2135 log << MSG::INFO << "in finalize(): Number of events in CgemClusterCreate" << m_totEvt << endreq;
2136 //delete myNTupleHelper;
2137
2138 return StatusCode::SUCCESS;
2139}
2140
2141void CgemClusterCreate::reset() {
2142 for (int i=0; i<3; i++) /* Layer: 3 */
2143 {
2144 for (int j=0; j<2; j++) /* Sheet: 2 */
2145 {
2146 for (int k=0; k<2; k++) /* Strip Type: 2 */
2147 {
2148 for (int l=0; l<1500; l++) /* StripID: 1419 */
2149 {
2150 m_strip[i][j][k][l] = -1;
2151 m_edep[i][j][k][l] = 0;
2152 }
2153 }
2154 }
2155 } /* End of 'for (int i=0; i<3; i++)' */
2156
2157 // reset the map
2158 m_x_map.clear();
2159 m_v_map.clear();
2160 m_xv_map.clear();
2161 m_trans_map.clear();
2162 m_driftrec_map.clear();
2163 m_pos_map.clear();
2164 m_strip_map.clear();
2165
2166}
2167
2168void CgemClusterCreate::resetFiredStripMap()
2169{
2170 for (int i=0; i<3; i++) /* Layer: 3 */
2171 {
2172 for (int j=0; j<2; j++) /* Sheet: 2 */
2173 {
2174 for (int k=0; k<2; k++) /* Strip Type: 2 */
2175 {
2176 myFiredStripMap[i][j][k].clear();
2177 }
2178 }
2179 }
2180
2181}
2182
2183/*
2184 void CgemClusterCreate::resetForCluster()
2185 {
2186 mySumQ=0.0;
2187 mySumQX=0.0;
2188 }
2189 */
2190
2191void CgemClusterCreate:: processstrips(int k){
2192
2193 for(int i=0;i<3;i++){
2194 for(int j=0;j<2;j++){
2195 ClusterFlag cluster;
2196 cluster.begin = 0;
2197 cluster.end = 0;
2198
2199 static int N_cluster;
2200 N_cluster = 0;
2201 for(int l=1;l<1500;l++){
2202 if(k==0){
2203 if(m_strip[i][j][k][l]>0&&m_strip[i][j][k][l-1]>0){
2204 // cout << "middle " << l << endl;
2205 cluster.end = l;}
2206 else{
2207 if(m_strip[i][j][k][l]<0&&m_strip[i][j][k][l-1]>0){
2208 // cout << "end" << l << endl;
2209 N_cluster= N_cluster + 1;
2210 // cout << "number" << N_cluster<< endl;
2211 m_x_group.push_back(cluster);
2212 //CgemCluster *cgemcluster = new CgemCluster();
2213 //cgemcluster->setlayerid(i);
2214 //cgemcluster->setsheetid(j);
2215 //cout << "sheet id" << cgemcluster->getsheetid() << endl;
2216 double energy=0;
2217 for(int n=cluster.begin;n<=cluster.end;n++){
2218 energy = energy + m_edep[i][j][k][n];
2219 }
2220 //cgemcluster->setenergydeposit(energy);
2221 //cgemcluster->setclusterid(N_cluster);
2222 //cgemcluster->setflag(k);
2223 //cgemcluster->setclusterflag(cluster.begin,cluster.end);
2224 RecCgemCluster *reccluster = new RecCgemCluster();
2225 reccluster->setlayerid(i);
2226 reccluster->setsheetid(j);
2227 // cout << "sheet id" << reccluster->getsheetid() << endl;
2228 reccluster->setenergydeposit(energy);
2229 reccluster->setclusterid(N_cluster);
2230 reccluster->setflag(k);
2231 reccluster->setclusterflag(cluster.begin,cluster.end);
2232 posxfind(reccluster);
2233 keytype key((10*i+j),reccluster->getclusterid());
2234 m_x_map[key] = reccluster;
2235 cluster.begin = l;
2236 cluster.end = l;
2237 }
2238 if(m_strip[i][j][k][l]>0&&m_strip[i][j][k][l-1]<0){
2239 // cout << "begin" << l << endl;
2240 cluster.begin = l;
2241 cluster.end = l;
2242 }
2243 }
2244 }
2245 if(k==1){
2246 if(m_strip[i][j][k][l]>0&&m_strip[i][j][k][l-1]>0){
2247 // cout << "middle " << l << endl;
2248 cluster.end = l;
2249 }
2250 else
2251 {
2252 if(m_strip[i][j][k][l]<0&&m_strip[i][j][k][l-1]>0){
2253 // cout << "end" << l << endl;
2254 N_cluster++;
2255 // cout << "number" << N_cluster<< endl;
2256 m_v_group.push_back(cluster);
2257 //CgemCluster* cgemcluster = new CgemCluster();
2258 //cgemcluster->setlayerid(i);
2259 //cgemcluster->setsheetid(j);
2260 double energy=0;
2261 for(int n=cluster.begin;n<=cluster.end;n++){
2262 energy = energy + m_edep[i][j][k][n];
2263 }
2264 //cgemcluster->setenergydeposit(energy);
2265 //cgemcluster->setclusterid(N_cluster);
2266 //cgemcluster->setflag(k);
2267 //cgemcluster->setclusterflag(cluster.begin,cluster.end);
2268 RecCgemCluster* reccluster = new RecCgemCluster();
2269 reccluster->setlayerid(i);
2270 reccluster->setsheetid(j);
2271 // cout << "sheet id" << reccluster->getsheetid() << endl;
2272 reccluster->setenergydeposit(energy);
2273 reccluster->setclusterid(N_cluster);
2274 reccluster->setflag(k);
2275 reccluster->setclusterflag(cluster.begin,cluster.end);
2276 posvfind(reccluster);
2277 keytype key((10*i+j),reccluster->getclusterid());
2278 m_v_map[key] = reccluster;
2279 cluster.begin = l;
2280 cluster.end = l;
2281 }
2282 if(m_strip[i][j][k][l]>0&&m_strip[i][j][k][l-1]<0){
2283 // cout << "begin" << l << endl;
2284 cluster.begin = l;
2285 cluster.end = l;
2286 }
2287 }
2288 }
2289 } //strip//
2290 m_x_group.clear();
2291 m_v_group.clear();
2292 } //sheet//
2293 } //layer//
2294} //end of the function//
2295
2296void CgemClusterCreate:: transcluster(){
2297 //convert v-cluster to x-direction//
2298 for(int i=0;i<3;i++){//layer 3//
2299 for(int j=0;j<2;j++){//sheet 2//
2300 for(int l=0;l<1500;l++){
2301
2302 keytype key((10*i+j),l);
2303 m_v_map_it = m_v_map.find(key);
2304 if(m_v_map_it != m_v_map.end()){
2305
2306 // cout << "found ID" << l << "in layer" << i << "in sheet" << j << endl;
2307 RecCgemCluster* reccluster = (*m_v_map_it).second;
2308 int b = reccluster->getclusterflagb();
2309 int e = reccluster->getclusterflage();
2310 // cout << "flag begin" << b << "flag end" << e << endl;
2311 ClusterFlag vcluster;
2312 vcluster.begin = b;
2313 vcluster.end = e;
2314 ClusterFlag cluster;
2315 double extrax = l_layer[i]*tan(a_stero[i]);
2316 cluster.begin = floor((b*W_pitch/(cos(a_stero[i]))-extrax)/W_pitch);
2317 cluster.end = floor(e/(cos(a_stero[i])));
2318 if(e*W_pitch/(cos(a_stero[i]))>w_sheet[i]){cluster.end = floor(w_sheet[i]/W_pitch);}
2319 if(cluster.begin<0){cluster.begin=0;}
2320 // cout << "cos angle" << cos(a_stero[i]) << "angle" << a_stero[i] <<endl;
2321 // cout << "trans begin" << cluster.begin <<endl;
2322 // cout << "trans end" << cluster.end <<endl;
2323 flagxv transflag;
2324 transflag.first = vcluster;
2325 transflag.second = cluster;
2326 m_trans_map[key]= transflag;
2327
2328 }
2329 }
2330 }
2331 }
2332}
2333
2334void CgemClusterCreate::mixcluster(){
2335
2336 for(int i=0;i<3;i++){
2337 for(int j=0;j<2;j++){
2338 for(int l=0;l<1500;l++){
2339
2340 keytype key((10*i+j),l);
2341 m_x_map_it = m_x_map.find(key);
2342 if(m_x_map_it != m_x_map.end()){
2343
2344
2345 // cout << "found x ID" << l << "in layer" << i << "in sheet" << j << endl;
2346 RecCgemCluster* xcluster = (*m_x_map_it).second;
2347 int xb = xcluster->getclusterflagb();
2348 int xe = xcluster->getclusterflage();
2349 // cout << "flag" << xb << xe << endl;
2350 for(int ll=0;ll<1500;ll++){
2351 keytype mkey((10*i+j),ll);
2352 m_trans_map_it = m_trans_map.find(mkey);
2353 if(m_trans_map_it == m_trans_map.end()){continue;}
2354 else{
2355 // cout << "found trans ID" << ll << "in layer" << i << "in sheet" << j << endl;
2356 flagxv transflag = (*m_trans_map_it).second;
2357 ClusterFlag vcluster = transflag.first;
2358 ClusterFlag cluster = transflag.second;
2359 if((xb>=cluster.begin&&xb<=cluster.end)||(cluster.begin>=xb&&cluster.begin<=xe)){
2360
2361 ClusterFlag x,v;
2362 x.begin = xb;
2363 x.end = xe;
2364 v.begin = cluster.begin;
2365 v.end = cluster.end;
2366 flagxv XV;
2367 XV.first = x;
2368 XV.second = v;
2369 m_xv_group.push_back(XV);
2370 m_mid_group.push_back(vcluster);
2371 RecCgemCluster* reccluster = new RecCgemCluster();
2372 reccluster->setlayerid(i);
2373 reccluster->setsheetid(j);
2374 reccluster->setflag(2);
2375 reccluster->setclusterid(m_xv_group.size());
2376 posxvfind(reccluster);
2377 keytype key((10*i+j),reccluster->getclusterid());
2378 m_xv_map[key] = reccluster;
2379 posindrift(reccluster);
2380 // cout << "###########" << endl;
2381 // cout << "xb" << xb << "xe" << xe << endl;
2382 // cout << "trans begin" << cluster.begin << "trans end" << cluster.end << endl;
2383 // cout << m_xv_group.size()<<endl;
2384 } //coincidence
2385
2386 else{continue;}
2387 } //find one trans cluster//
2388 } //loop m_trans_map//
2389 }// find x cluster//
2390 }//cluster//
2391 }//sheet//
2392 }//layer//
2393 m_xv_group.clear();
2394 m_mid_group.clear();
2395}//end of function mixcluster//
2396
2397
2398void CgemClusterCreate::posxfind(RecCgemCluster* reccluster){
2399
2400 double phi = 0.0;
2401 int layerid = reccluster->getlayerid();
2402 int sheetid = reccluster->getsheetid();
2403 int begin = reccluster->getclusterflagb();
2404 int end = reccluster->getclusterflage();
2405 for(int k=begin;k<=end;k++){
2406 phi = phi + m_edep[layerid][sheetid][0][k]*k*W_pitch;
2407 }
2408 reccluster->setrecphi(phi/(reccluster->getenergydeposit()));
2409}
2410
2411void CgemClusterCreate::posvfind(RecCgemCluster* reccluster){
2412 double v = 0.0;
2413 int layerid = reccluster->getlayerid();
2414 int sheetid = reccluster->getsheetid();
2415 int begin = reccluster->getclusterflagb();
2416 int end = reccluster->getclusterflage();
2417 for(int k=begin;k<=end;k++){
2418 v = v + m_edep[layerid][sheetid][1][k]*k*W_pitch;
2419 }
2420 reccluster->setrecv(v/(reccluster->getenergydeposit()));
2421}
2422
2423void CgemClusterCreate::posxvfind(RecCgemCluster* reccluster){
2424 int layerid = reccluster->getlayerid();
2425 int sheetid = reccluster->getsheetid();
2426 int clusterid = reccluster->getclusterid();
2427 ClusterFlag x = (m_xv_group[clusterid-1]).first;;
2428 ClusterFlag v = (m_xv_group[clusterid-1]).second;
2429 ClusterFlag mid = m_mid_group[clusterid-1];
2430 //cout << "&&&&&&x" << x.begin << x.end << endl;
2431 //cout << "&&&&&&v" << v.begin << v.end << endl;
2432 //cout << "&&&&&&mid" << mid.begin << mid.end << endl;
2433 double phi = 0.0;
2434 double recv = 0.0;
2435 double xenergy = 0.0;
2436 double venergy = 0.0;
2437 double pphi = 0.0;
2438 double pv = 0.0;
2439 for(int ii=x.begin;ii<=x.end;ii++){
2440 for(int jj=v.begin;jj<=v.end;jj++){
2441 if(ii==jj){
2442 // cout << "same ID" << ii << jj << endl;
2443 xenergy = xenergy + m_edep[layerid][sheetid][0][ii];
2444 phi = phi + m_edep[layerid][sheetid][0][ii]*ii*W_pitch;
2445 pphi= pphi+ii*W_pitch;
2446 m_sameID.push_back(ii);
2447 //cout << "!!!!!!!!!!1energy" << xenergy << endl;
2448 //cout << "!!!!!!!!!!phi" << phi << endl;
2449 }
2450 else{continue;}
2451 }
2452 }
2453 int nvstrip=0;
2454 for(int kk=mid.begin;kk<=mid.end;kk++){
2455 double extrax = l_layer[layerid]*tan(a_stero[layerid]);
2456 double pa = (kk*W_pitch/(cos(a_stero[layerid]))-extrax)/W_pitch;
2457 double pb = kk/(cos(a_stero[layerid]));
2458 if(pb>(w_sheet[layerid]/W_pitch)){pb = w_sheet[layerid]/W_pitch;}
2459 if(pa<0){pa = 0;}
2460 for(int ll=0;ll<(int)m_sameID.size();ll++){
2461 if(pa<=m_sameID[ll]&&pb>=m_sameID[ll]){
2462 // cout << "pa" << pa << "pb" << pb << endl;
2463 venergy = venergy + m_edep[layerid][sheetid][1][kk];
2464 recv = recv + m_edep[layerid][sheetid][1][kk]*kk*W_pitch;
2465 pv = pv + kk*W_pitch;
2466 nvstrip++;
2467 //cout << "$$$$$$$$$pv" << pv << endl;
2468 //cout << "!!!!!!!!!kk" << kk << endl;
2469 //cout << "!!!!!!!!!2energy" << venergy << endl;
2470 //cout << "!!!!!!!!!recv" << recv << endl;
2471 break;
2472 }
2473 else{continue;}
2474 }
2475 }
2476 int nxstrip = m_sameID.size();
2477 //cout << "!!!!!!!!pvall " << pv << endl;
2478 //cout << "!!!!!!!!xstrip" << nxstrip << endl;
2479 //cout << "!!!!!!!!vstrip" << nvstrip << endl;
2480 keytype key((10*layerid+sheetid),clusterid);
2481 pphi = pphi/nxstrip;
2482 pv = pv/nvstrip;
2483 keytype strip(nxstrip,nvstrip);
2484 postype pos(pphi,pv);
2485 //cout << "!!!!!!!!pphi" << pphi << endl;
2486 //cout << "!!!!!!!!pv" << pv << endl;
2487 m_strip_map[key]=strip;
2488 m_pos_map[key]=pos;
2489 m_sameID.clear();
2490 reccluster->setenergydeposit(xenergy+venergy);
2491 //cout << "!!!!get phi" << phi/xenergy << endl;
2492 //cout << "!!!!get v " << recv/venergy << endl;
2493 reccluster->setrecphi(phi/xenergy);
2494 reccluster->setrecv(recv/venergy);
2495}
2496
2497void CgemClusterCreate::posindrift(RecCgemCluster* reccluster){
2498 int layerid = reccluster->getlayerid();
2499 int sheetid = reccluster->getsheetid();
2500 int clusterid = reccluster->getclusterid();
2501 double dphi = ((reccluster->getrecphi())+sheetid*w_sheet[layerid])/r_layer[layerid];
2502 double dz = ((reccluster->getrecv())-(reccluster->getrecphi())*cos(a_stero[layerid]))/sin(a_stero[layerid]);
2503 double posphi = dr_layer[layerid]*dphi;
2504 for(int vv =0;vv<n_sheet[layerid];vv++){
2505 if (posphi>= vv*dw_sheet[layerid] && posphi < (vv+1)*dw_sheet[layerid]){
2506 posphi = posphi-vv*dw_sheet[layerid];
2507 }
2508 }
2509 double posv = posphi*cos(a_stero[layerid])+dz*sin(a_stero[layerid]);
2510 //cout << "!!!!!!!!dphi" << dphi << endl;
2511 //cout << "!!!!!!!!dr" << dr_layer[layerid]*dphi << endl;
2512 //cout << "!!!!!!!!posphi" << posphi << endl;
2513 //cout << "!!!!!!!!posv " << posv << endl;
2514 postype positiond(posphi,posv);
2515 keytype key((10*layerid+sheetid),clusterid);
2516 m_driftrec_map[key] = positiond;
2517}
2518
2519void CgemClusterCreate::fillMCTruth()
2520{
2521 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
2522 if (!mcParticleCol)
2523 {
2524 std::cout << "Could not retrieve McParticelCol" << std::endl;
2525 return;
2526 }
2527 int i_mc=0;
2528 int foundMom=0;
2529 int startDecay=0;
2530 int itmp=0;
2531 int mcPDG[100];
2532 double Pmc[100],costhmc[100],phimc[100];
2533 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
2534 for (; iter_mc != mcParticleCol->end(); iter_mc++)
2535 {
2536 int idx = (*iter_mc)->trackIndex();
2537 int pdgid = (*iter_mc)->particleProperty();
2538 int mother_pdg = ((*iter_mc)->mother()).particleProperty();
2539 int mother_idx = (*iter_mc)->mother().trackIndex();
2540 int primaryParticle = (*iter_mc)->primaryParticle();
2541 int fromGenerator = (*iter_mc)->decayFromGenerator();
2542 if(primaryParticle==1) continue;
2543 if(fromGenerator==0) continue;
2544 if(pdgid==myMotherParticleID) {
2545 foundMom=idx;
2546 continue;
2547 }
2548 if(foundMom==0) continue;
2549 if(mother_pdg==myMotherParticleID) {
2550 startDecay=idx;
2551
2552 }
2553 if(pdgid!=myMotherParticleID&&foundMom>0&&startDecay==0)
2554 {
2555 itmp++;
2556 continue;
2557 }
2558 mother_idx=mother_idx-foundMom-itmp;
2559
2560 if(myPrintFlag==1) {
2561 if(i_mc==0) {
2562 cout<<"****** MC particles ****** "<<endl;
2563 cout <<setw(10)<<"index"
2564 <<setw(10)<<"pdg"
2565 <<setw(16)<<"primaryParticle"
2566 <<setw(15)<<"fromGenerator"
2567 <<setw(15)<<"from_trk"
2568 <<endl;
2569 }
2570 cout <<setw(10)<<i_mc
2571 <<setw(10)<<pdgid
2572 <<setw(16)<<primaryParticle
2573 <<setw(15)<<fromGenerator
2574 <<setw(15)<<mother_idx
2575 <<endl;
2576 }
2577 if(i_mc<100)
2578 {
2579 mcPDG[i_mc]=pdgid;
2580 HepLorentzVector p4_mc = (*iter_mc)->initialFourMomentum();
2581 Pmc[i_mc]=p4_mc.rho();
2582 costhmc[i_mc]=p4_mc.cosTheta();
2583 phimc[i_mc]=p4_mc.phi();
2584 }
2585 i_mc++;
2586 }// loop MC particles
2587 if(i_mc>=100) i_mc=100;
2588
2589 int nTruth[3]={0,0,0};
2590 int trkID_Layer1[100], trkID_Layer2[100], trkID_Layer3[100];
2591 int pdg_Layer1[100], pdg_Layer2[100], pdg_Layer3[100];
2592 double x_pre_Layer1[100], x_pre_Layer2[100], x_pre_Layer3[100];
2593 double y_pre_Layer1[100], y_pre_Layer2[100], y_pre_Layer3[100];
2594 double z_pre_Layer1[100], z_pre_Layer2[100], z_pre_Layer3[100];
2595 double x_post_Layer1[100], x_post_Layer2[100], x_post_Layer3[100];
2596 double y_post_Layer1[100], y_post_Layer2[100], y_post_Layer3[100];
2597 double z_post_Layer1[100], z_post_Layer2[100], z_post_Layer3[100];
2598 double phi_Layer1[100], phi_Layer2[100], phi_Layer3[100];
2599 double z_Layer1[100], z_Layer2[100], z_Layer3[100];
2600 double V_Layer1[100], V_Layer2[100], V_Layer3[100];
2601 double px_pre_Layer1[100], px_pre_Layer2[100], px_pre_Layer3[100];
2602 double py_pre_Layer1[100], py_pre_Layer2[100], py_pre_Layer3[100];
2603 double pz_pre_Layer1[100], pz_pre_Layer2[100], pz_pre_Layer3[100];
2604 double en_Layer1[100], en_Layer2[100], en_Layer3[100];
2605 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
2606 if (!cgemMcHitCol)
2607 {
2608 std::cout << "Could not retrieve cgemMcHitCol" << std::endl;
2609 return;
2610 }
2611 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
2612 for (; iter_truth!= cgemMcHitCol->end(); iter_truth++)
2613 {
2614 int iLayer = (*iter_truth)->GetLayerID();
2615 int trkID = (*iter_truth)->GetTrackID();
2616 int pdg = (*iter_truth)->GetPDGCode();
2617 double x_pre = (*iter_truth)->GetPositionXOfPrePoint();// in mm
2618 double y_pre = (*iter_truth)->GetPositionYOfPrePoint();
2619 double z_pre = (*iter_truth)->GetPositionZOfPrePoint();
2620 double x_post = (*iter_truth)->GetPositionXOfPostPoint();
2621 double y_post = (*iter_truth)->GetPositionYOfPostPoint();
2622 double z_post = (*iter_truth)->GetPositionZOfPostPoint();
2623 double x_middle = 0.5*(x_pre+x_post);
2624 double y_middle = 0.5*(y_pre+y_post);
2625 double z_middle = 0.5*(z_pre+z_post);
2626 double r_middle = sqrt(x_middle*x_middle+y_middle*y_middle); //cout<<"r_middle = "<<r_middle<<endl;
2627 double phi_middle = atan2(y_middle, x_middle);
2628 double px_pre = (*iter_truth)->GetMomentumXOfPrePoint();
2629 double py_pre = (*iter_truth)->GetMomentumYOfPrePoint();
2630 double pz_pre = (*iter_truth)->GetMomentumZOfPrePoint();
2631 double en = (*iter_truth)->GetTotalEnergyDeposit();
2632 CgemGeoReadoutPlane* readoutPlane=NULL;
2633 //cout<<"phi_middle, z_middle = "<<phi_middle<<", "<<z_middle<<endl;
2634 //cout<<"iLayer = "<<iLayer<<endl;
2635 for(int j=0; j<myNSheets[iLayer]; j++)
2636 {
2637 //cout<<"j = "<<j <<endl;
2638 readoutPlane=myCgemGeomSvc->getReadoutPlane(iLayer,j);
2639 //cout<<"OnThePlane: "<<readoutPlane->OnThePlane(phi_middle,z_middle)<<endl;
2640 if(readoutPlane->OnThePlane(phi_middle,z_middle)) break;
2641 readoutPlane=NULL;
2642 }
2643 double V_mid = 9999;
2644 if(readoutPlane==NULL) {
2645 if(myPrintFlag) cout<<"CgemClusterCreate::fillMCTruth() readoutPlane not found with phi_middle = "<<phi_middle<<", layer = "<<iLayer<<endl;
2646 }
2647 else {
2648 //readoutPlane->print();
2649 //cout<<"phi_middle, z_middle = "<<phi_middle<<", "<<z_middle<<endl;
2650 V_mid = readoutPlane->getVFromPhiZ(phi_middle,z_middle);
2651 //cout<<"V_mid = "<<V_mid<<endl;
2652 //double V_mid = readoutPlane->GetVFromLocalXZ(phi_middle,z_middle);
2653 }
2654
2655 switch(iLayer)
2656 {
2657 case 0:
2658 if(nTruth[0]>=100) break;
2659 trkID_Layer1[nTruth[0]]=trkID;
2660 pdg_Layer1[nTruth[0]]=pdg;
2661 x_pre_Layer1[nTruth[0]]=x_pre;
2662 y_pre_Layer1[nTruth[0]]=y_pre;
2663 z_pre_Layer1[nTruth[0]]=z_pre;
2664 x_post_Layer1[nTruth[0]]=x_post;
2665 y_post_Layer1[nTruth[0]]=y_post;
2666 z_post_Layer1[nTruth[0]]=z_post;
2667 phi_Layer1[nTruth[0]]=atan2(y_post+y_pre,x_post+x_pre);
2668 z_Layer1[nTruth[0]]=z_middle;
2669 V_Layer1[nTruth[0]]=V_mid;
2670 px_pre_Layer1[nTruth[0]]=px_pre;
2671 py_pre_Layer1[nTruth[0]]=py_pre;
2672 pz_pre_Layer1[nTruth[0]]=pz_pre;
2673 en_Layer1[nTruth[0]]=en;
2674 break;
2675 case 1:
2676 if(nTruth[1]>=100) break;
2677 trkID_Layer2[nTruth[1]]=trkID;
2678 pdg_Layer2[nTruth[1]]=pdg;
2679 x_pre_Layer2[nTruth[1]]=x_pre;
2680 y_pre_Layer2[nTruth[1]]=y_pre;
2681 z_pre_Layer2[nTruth[1]]=z_pre;
2682 x_post_Layer2[nTruth[1]]=x_post;
2683 y_post_Layer2[nTruth[1]]=y_post;
2684 z_post_Layer2[nTruth[1]]=z_post;
2685 phi_Layer2[nTruth[1]]=atan2(y_post+y_pre,x_post+x_pre);
2686 z_Layer2[nTruth[1]]=z_middle;
2687 V_Layer2[nTruth[1]]=V_mid;
2688 px_pre_Layer2[nTruth[1]]=px_pre;
2689 py_pre_Layer2[nTruth[1]]=py_pre;
2690 pz_pre_Layer2[nTruth[1]]=pz_pre;
2691 en_Layer2[nTruth[1]]=en;
2692 break;
2693 case 2:
2694 if(nTruth[2]>=100) break;
2695 trkID_Layer3[nTruth[2]]=trkID;
2696 pdg_Layer3[nTruth[2]]=pdg;
2697 x_pre_Layer3[nTruth[2]]=x_pre;
2698 y_pre_Layer3[nTruth[2]]=y_pre;
2699 z_pre_Layer3[nTruth[2]]=z_pre;
2700 x_post_Layer3[nTruth[2]]=x_post;
2701 y_post_Layer3[nTruth[2]]=y_post;
2702 z_post_Layer3[nTruth[2]]=z_post;
2703 phi_Layer3[nTruth[2]]=atan2(y_post+y_pre,x_post+x_pre);
2704 z_Layer3[nTruth[2]]=z_middle;
2705 V_Layer3[nTruth[2]]=V_mid;
2706 px_pre_Layer3[nTruth[2]]=px_pre;
2707 py_pre_Layer3[nTruth[2]]=py_pre;
2708 pz_pre_Layer3[nTruth[2]]=pz_pre;
2709 en_Layer3[nTruth[2]]=en;
2710 break;
2711 default:
2712 cout<<"wrong layer ID for CGEM: "<<iLayer<<endl;
2713 }
2714 nTruth[iLayer]++;
2715 }// loop hit truth
2716 for(int i=0; i<3; i++) if(nTruth[i]>100) nTruth[i]=100;
2717
2718 if(myNtuple)
2719 {
2720 myNTupleHelper->fillArrayInt("pdgmc","nmc",(int*)mcPDG,i_mc);
2721 myNTupleHelper->fillArray("pmc","nmc",(double*)Pmc,i_mc);
2722 myNTupleHelper->fillArray("costhmc","nmc",(double*)costhmc,i_mc);
2723 myNTupleHelper->fillArray("phimc","nmc",(double*)phimc,i_mc);
2724
2725 myNTupleHelper->fillArrayInt("trkID_Layer1","nTruthLay1",(int*) trkID_Layer1, nTruth[0]);
2726 myNTupleHelper->fillArrayInt( "pdg_Layer1","nTruthLay1",(int*) pdg_Layer1, nTruth[0]);
2727 myNTupleHelper->fillArray( "x_pre_Layer1","nTruthLay1",(double*) x_pre_Layer1, nTruth[0]);
2728 myNTupleHelper->fillArray( "y_pre_Layer1","nTruthLay1",(double*) y_pre_Layer1, nTruth[0]);
2729 myNTupleHelper->fillArray( "z_pre_Layer1","nTruthLay1",(double*) z_pre_Layer1, nTruth[0]);
2730 myNTupleHelper->fillArray( "x_post_Layer1","nTruthLay1",(double*) x_post_Layer1, nTruth[0]);
2731 myNTupleHelper->fillArray( "y_post_Layer1","nTruthLay1",(double*) y_post_Layer1, nTruth[0]);
2732 myNTupleHelper->fillArray( "z_post_Layer1","nTruthLay1",(double*) z_post_Layer1, nTruth[0]);
2733 myNTupleHelper->fillArray( "px_pre_Layer1","nTruthLay1",(double*) px_pre_Layer1, nTruth[0]);
2734 myNTupleHelper->fillArray( "py_pre_Layer1","nTruthLay1",(double*) py_pre_Layer1, nTruth[0]);
2735 myNTupleHelper->fillArray( "pz_pre_Layer1","nTruthLay1",(double*) pz_pre_Layer1, nTruth[0]);
2736 myNTupleHelper->fillArray( "en_Layer1","nTruthLay1",(double*) en_Layer1, nTruth[0]);
2737 myNTupleHelper->fillArray( "phi_Layer1","nTruthLay1",(double*) phi_Layer1, nTruth[0]);
2738 myNTupleHelper->fillArray( "V_Layer1","nTruthLay1",(double*) V_Layer1, nTruth[0]);
2739
2740 myNTupleHelper->fillArrayInt("trkID_Layer2","nTruthLay2",(int*) trkID_Layer2, nTruth[1]);
2741 myNTupleHelper->fillArrayInt( "pdg_Layer2","nTruthLay2",(int*) pdg_Layer2, nTruth[1]);
2742 myNTupleHelper->fillArray( "x_pre_Layer2","nTruthLay2",(double*) x_pre_Layer2, nTruth[1]);
2743 myNTupleHelper->fillArray( "y_pre_Layer2","nTruthLay2",(double*) y_pre_Layer2, nTruth[1]);
2744 myNTupleHelper->fillArray( "z_pre_Layer2","nTruthLay2",(double*) z_pre_Layer2, nTruth[1]);
2745 myNTupleHelper->fillArray( "x_post_Layer2","nTruthLay2",(double*) x_post_Layer2, nTruth[1]);
2746 myNTupleHelper->fillArray( "y_post_Layer2","nTruthLay2",(double*) y_post_Layer2, nTruth[1]);
2747 myNTupleHelper->fillArray( "z_post_Layer2","nTruthLay2",(double*) z_post_Layer2, nTruth[1]);
2748 myNTupleHelper->fillArray( "px_pre_Layer2","nTruthLay2",(double*) px_pre_Layer2, nTruth[1]);
2749 myNTupleHelper->fillArray( "py_pre_Layer2","nTruthLay2",(double*) py_pre_Layer2, nTruth[1]);
2750 myNTupleHelper->fillArray( "pz_pre_Layer2","nTruthLay2",(double*) pz_pre_Layer2, nTruth[1]);
2751 myNTupleHelper->fillArray( "en_Layer2","nTruthLay2",(double*) en_Layer2, nTruth[1]);
2752 myNTupleHelper->fillArray( "phi_Layer2","nTruthLay2",(double*) phi_Layer2, nTruth[1]);
2753 myNTupleHelper->fillArray( "V_Layer2","nTruthLay2",(double*) V_Layer2, nTruth[1]);
2754
2755 myNTupleHelper->fillArrayInt("trkID_Layer3","nTruthLay3",(int*) trkID_Layer3, nTruth[2]);
2756 myNTupleHelper->fillArrayInt( "pdg_Layer3","nTruthLay3",(int*) pdg_Layer3, nTruth[2]);
2757 myNTupleHelper->fillArray( "x_pre_Layer3","nTruthLay3",(double*) x_pre_Layer3, nTruth[2]);
2758 myNTupleHelper->fillArray( "y_pre_Layer3","nTruthLay3",(double*) y_pre_Layer3, nTruth[2]);
2759 myNTupleHelper->fillArray( "z_pre_Layer3","nTruthLay3",(double*) z_pre_Layer3, nTruth[2]);
2760 myNTupleHelper->fillArray( "x_post_Layer3","nTruthLay3",(double*) x_post_Layer3, nTruth[2]);
2761 myNTupleHelper->fillArray( "y_post_Layer3","nTruthLay3",(double*) y_post_Layer3, nTruth[2]);
2762 myNTupleHelper->fillArray( "z_post_Layer3","nTruthLay3",(double*) z_post_Layer3, nTruth[2]);
2763 myNTupleHelper->fillArray( "px_pre_Layer3","nTruthLay3",(double*) px_pre_Layer3, nTruth[2]);
2764 myNTupleHelper->fillArray( "py_pre_Layer3","nTruthLay3",(double*) py_pre_Layer3, nTruth[2]);
2765 myNTupleHelper->fillArray( "pz_pre_Layer3","nTruthLay3",(double*) pz_pre_Layer3, nTruth[2]);
2766 myNTupleHelper->fillArray( "en_Layer3","nTruthLay3",(double*) en_Layer3, nTruth[2]);
2767 myNTupleHelper->fillArray( "phi_Layer3","nTruthLay3",(double*) phi_Layer3, nTruth[2]);
2768 myNTupleHelper->fillArray( "V_Layer3","nTruthLay3",(double*) V_Layer3, nTruth[2]);
2769 }
2770
2771}
2772
2773void CgemClusterCreate::checkRecCgemClusterCol()
2774{
2775 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),"/Event/Recon/RecCgemClusterCol");
2776 if(aCgemClusterCol)
2777 {
2778 RecCgemClusterCol::iterator iter_cluster=aCgemClusterCol->begin();
2779 int nCluster = aCgemClusterCol->size();
2780 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~ check RecCgemClusterCol:"<<endl;
2781 cout <<setw(10)<<"idx"
2782 <<setw(10)<<"layer"
2783 <<setw(10)<<"sheet"
2784 <<setw(10)<<"XVFlag"
2785 <<setw(10)<<"id1 ~"
2786 <<setw(10)<<"id2"
2787 <<setw(15)<<"phi"
2788 <<setw(15)<<"V"
2789 <<setw(15)<<"z"
2790 <<endl;
2791 for(; iter_cluster!=aCgemClusterCol->end(); iter_cluster++)
2792 {
2793 cout <<setw(10)<<(*iter_cluster)->getclusterid()
2794 <<setw(10)<<(*iter_cluster)->getlayerid()
2795 <<setw(10)<<(*iter_cluster)->getsheetid()
2796 <<setw(10)<<(*iter_cluster)->getflag()
2797 <<setw(10)<<(*iter_cluster)->getclusterflagb()
2798 <<setw(10)<<(*iter_cluster)->getclusterflage()
2799 <<setw(15)<<setprecision(10)<<(*iter_cluster)->getrecphi()
2800 <<setw(15)<<setprecision(10)<<(*iter_cluster)->getrecv()
2801 <<setw(15)<<setprecision(10)<<(*iter_cluster)->getRecZ()
2802 <<endl;
2803 }
2804 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
2805 }
2806 else cout<<"CgemClusterCreate::checkRecCgemClusterCol(): RecCgemClusterCol not accessible!"<<endl;
2807}
2808
2809bool CgemClusterCreate::isGoodDigi(CgemDigiCol::iterator iter)
2810{
2811 double time = (*iter)->getTime_ns();
2812 bool isGood=true;
2813 if(time<-180.||time>100.) return false;
2814 double Q = (*iter)->getCharge_fc();
2815 if(Q<0.) return false;
2816 return isGood;
2817}
2818
2819
2820bool CgemClusterCreate::selDigi(CgemDigiCol::iterator iter, int sel)
2821{
2822 if(sel==0) return true;
2823 else {
2824 double time = (*iter)->getTime_ns();
2825 bool timeIsGood=true;
2826 if(time<m_minDigiTime||time>m_maxDigiTime) timeIsGood=false;
2827
2828 double Q = (*iter)->getCharge_fc();
2829 bool chargeIsGood=true;
2830 if(Q<myQMin) chargeIsGood=false;
2831
2832 const Identifier ident = (*iter)->identify();
2833 int layer = CgemID::layer(ident);
2834 int sheet = CgemID::sheet(ident);
2835 int strip = CgemID::strip(ident);
2836 int view = 1;
2837 bool is_xstrip = CgemID::is_xstrip(ident);
2838 if(is_xstrip == true) view = 0;
2839 int quality = lutreader->GetQuality(layer, sheet, view, strip);
2840 bool qualityIsGood=true;
2841 if(quality!=0) qualityIsGood=false;
2842
2843 if(sel==1) return timeIsGood&&chargeIsGood;
2844 else if(sel==-1) return !timeIsGood&&chargeIsGood;
2845 else if(sel==2) return timeIsGood&&chargeIsGood&&qualityIsGood;
2846 }
2847 return false;
2848
2849}
2850
2851float CgemClusterCreate::get_Time(CgemDigiCol::iterator iDigiCol){
2852 //Get digi time
2853 float time = (*iDigiCol)->getTime_ns();
2854 //Get rising time from calibration
2855 float time_rising = get_TimeRising(iDigiCol);
2856 //Get time-walk
2857 float time_walk = get_TimeWalk(iDigiCol);
2858 //Get time-reference
2859 float time_reference = get_TimeReference(iDigiCol);
2860 //
2861 float time_shift_custom = -35;
2862 //cout<<time<<" "<<time_rising<<" "<<time_walk<<" "<<time_reference;
2863 time-=(time_rising+time_walk+time_reference+time_shift_custom);
2864 //cout<<" "<<time<<endl;
2865 return time;
2866}
2867
2868float CgemClusterCreate::get_TimeRising(CgemDigiCol::iterator iDigiCol){
2869 float time_rising=0;
2870 //Get the digi information
2871 const Identifier ident = (*iDigiCol)->identify();
2872 int layerid = CgemID::layer(ident);
2873 int sheetid = CgemID::sheet(ident);
2874 int stripid = CgemID::strip(ident);
2875 int view = 1;
2876 bool is_xstrip = CgemID::is_xstrip(ident);
2877 if(is_xstrip == true) view = 0;
2878 float charge = (*iDigiCol)->getCharge_fc();
2879 //Get rising time from calibration
2880 time_rising = myCgemCalibSvc->getTimeRising(layerid, view, sheetid, stripid, charge, 0.);
2881 return time_rising;
2882}
2883
2884float CgemClusterCreate::get_TimeWalk(CgemDigiCol::iterator iDigiCol){
2885 float time_walk=0;
2886 const Identifier ident = (*iDigiCol)->identify();
2887 int layerid = CgemID::layer(ident);
2888 int sheetid = CgemID::sheet(ident);
2889 int stripid = CgemID::strip(ident);
2890 int view = 1;
2891 bool is_xstrip = CgemID::is_xstrip(ident);
2892 if(is_xstrip == true) view = 0;
2893 float thr = lutreader->Get_thr_T_fC(layerid, sheetid, view, stripid);
2894 float charge = (*iDigiCol)->getChargeChannel();
2895 time_walk = myCgemCalibSvc->getTimeWalk(charge,thr);
2896 return time_walk;
2897}
2898
2899float CgemClusterCreate::get_TimeReference(CgemDigiCol::iterator iDigiCol){
2900 float time_reference=0;
2901 const Identifier ident = (*iDigiCol)->identify();
2902 int layerid = CgemID::layer(ident);
2903 int sheetid = CgemID::sheet(ident);
2904 int stripid = CgemID::strip(ident);
2905 int view = 1;
2906 bool is_xstrip = CgemID::is_xstrip(ident);
2907 if(is_xstrip == true) view = 0;
2908 time_reference += lutreader->GetSignal_FEBStartTime_ns(layerid, sheetid, view, stripid);
2909 time_reference += lutreader->GetSignal_StartTime_ns(layerid, sheetid, view, stripid);
2910 return time_reference;
2911}
tan
double tan(const BesAngle a)
Definition BesAngle.h:216
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
dr_layer
const double dr_layer[3]
Definition CgemClusterCreate.cxx:56
drift_gap
const double drift_gap
Definition CgemClusterCreate.cxx:70
n_sheet
const int n_sheet[3]
Definition CgemClusterCreate.cxx:58
a_stero
const double a_stero[3]
Definition CgemClusterCreate.cxx:52
dw_sheet
const double dw_sheet[3]
Definition CgemClusterCreate.cxx:54
r_layer
const double r_layer[3]
Definition CgemClusterCreate.cxx:55
W_pitch
#define W_pitch
Definition CgemClusterCreate.cxx:46
l_layer
const int l_layer[3]
Definition CgemClusterCreate.cxx:57
w_sheet
const double w_sheet[3]
Definition CgemClusterCreate.cxx:53
n
const Int_t n
Definition DataBase/tau_mode.c:65
phi2
Double_t phi2
Definition DataBase/tau_mode.c:7
f1
TFile * f1
Definition DataBase/tau_mode.c:4
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
phi1
Double_t phi1
Definition DataBase/tau_mode.c:7
time
Double_t time
Definition DataBase/tau_mode.c:7
abs
double abs(const EvtComplex &c)
Definition EvtComplex.hh:212
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_nphot INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition KK2f.h:50
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
RecCgemClusterCol
ObjectVector< RecCgemCluster > RecCgemClusterCol
Definition RecCgemCluster.h:104
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
M_PI
#define M_PI
Definition TConstant.h:4
key
*************DOUBLE PRECISION m_pi *DOUBLE PRECISION m_HvecTau2 DOUBLE PRECISION m_HvClone2 DOUBLE PRECISION m_gamma1 DOUBLE PRECISION m_gamma2 DOUBLE PRECISION m_thet1 DOUBLE PRECISION m_thet2 INTEGER m_IFPHOT *COMMON c_Taupair $ !Spin Polarimeter vector first Tau $ !Spin Polarimeter vector second Tau $ !Clone Spin Polarimeter vector first Tau $ !Clone Spin Polarimeter vector second Tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !phi of HvecTau1 $ !theta of HvecTau1 $ !phi of HvecTau2 $ !theta of HvecTau2 $ !super key
Definition Taupair.h:42
BesTimer::start
void start(void)
Definition BesTimer.cxx:27
BesTimer::elapsed
float elapsed(void) const
Definition BesTimer.h:23
BesTimer::stop
void stop(void)
Definition BesTimer.cxx:39
CgemClusterCreate::CgemClusterCreate
CgemClusterCreate(const std::string &name, ISvcLocator *pSvcLocator)
Definition CgemClusterCreate.cxx:73
CgemGeoLayer::getOuterROfAnodeCu1
double getOuterROfAnodeCu1() const
Definition CgemGeoLayer.h:164
CgemGeoLayer::getInnerROfAnodeCu1
double getInnerROfAnodeCu1() const
Definition CgemGeoLayer.h:163
CgemGeoLayer::getNumberOfSheet
int getNumberOfSheet() const
Definition CgemGeoLayer.h:215
CgemGeoLayer::getNumberOfChannelX
int getNumberOfChannelX() const
Definition CgemGeoLayer.h:216
CgemGeoLayer::getOrientation
bool getOrientation() const
Definition CgemGeoLayer.h:27
CgemGeoLayer::getOuterROfAnodeCu2
double getOuterROfAnodeCu2() const
Definition CgemGeoLayer.h:168
CgemGeoLayer::getInnerROfAnodeCu2
double getInnerROfAnodeCu2() const
Definition CgemGeoLayer.h:167
CgemGeoReadoutPlane::getVFromPhiZ
double getVFromPhiZ(double phi, double z, bool checkRange=true) const
Definition CgemGeoReadoutPlane.h:178
CgemGeoReadoutPlane::getVIDInNextSheetFromVID
int getVIDInNextSheetFromVID(int vID, double phimin_next) const
Definition CgemGeoReadoutPlane.h:200
CgemGeoReadoutPlane::getZFromPhiV
double getZFromPhiV(double phi, double V, int checkXRange=1) const
Definition CgemGeoReadoutPlane.h:213
CgemGeoReadoutPlane::getCentralVFromVID
double getCentralVFromVID(int V_ID) const
Definition CgemGeoReadoutPlane.h:78
CgemGeoReadoutPlane::getPhiFromXID
double getPhiFromXID(int X_ID) const
Definition CgemGeoReadoutPlane.cxx:267
CgemGeoReadoutPlane::OnThePlane
bool OnThePlane(double phi, double z) const
Definition CgemGeoReadoutPlane.h:154
CgemGeoReadoutPlane::getVInNextSheetFromV
double getVInNextSheetFromV(double v, double phiminNext) const
Definition CgemGeoReadoutPlane.h:61
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition CgemGeomSvc.h:48
CgemGeomSvc::getNumberOfCgemLayer
int getNumberOfCgemLayer() const
Definition CgemGeomSvc.h:45
CgemGeomSvc::getReadoutPlane
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition CgemGeomSvc.h:140
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::GetSignal_StartTime_ns
float GetSignal_StartTime_ns(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:342
CgemLUTReader::GetQuality
int GetQuality(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:378
CgemLUTReader::Get_thr_T_fC
float Get_thr_T_fC(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:282
CgemLUTReader::GetSignal_FEBStartTime_ns
float GetSignal_FEBStartTime_ns(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:360
IBesTimerSvc::addItem
virtual BesTimer * addItem(const std::string &name)=0
ICgemCalibFunSvc::getSigma
virtual double getSigma(int layer, int xvFlag, int readoutMode, double angle, double Q, double T) const =0
ICgemCalibFunSvc::getTimeRising
virtual double getTimeRising(int layer, int xvFlag, int sheet, int stripID, double Q=100., double z=0.) const =0
ICgemCalibFunSvc::getTimeFalling
virtual double getTimeFalling(int layer, int xvFlag, int sheet, int stripID, double Q=100., double z=0.) const =0
ICgemCalibFunSvc::getTimeWalk
virtual double getTimeWalk(int layer, int xvFlag, int sheet, int stripID, double Q) const =0
ICgemGeomSvc::getReadoutPlane
virtual CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const =0
KalmanFit::Helix::dr
double dr(void) const
returns an element of parameters.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:229
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Reconstruction/KalFitAlg/KalFitAlg-00-15-20/KalFitAlg/helix/Helix.h:211
NTupleHelper::fillLong
void fillLong(string name, long value)
Definition NTupleHelper.cxx:86
NTupleHelper::fillDouble
void fillDouble(string name, double value)
Definition NTupleHelper.cxx:75
NTupleHelper::fillArrayInt
void fillArrayInt(string name, string index_name, int *value, int size)
Definition NTupleHelper.cxx:157
NTupleHelper::fillArray
void fillArray(string name, string index_name, double *value, int size)
Definition NTupleHelper.cxx:99
RecCgemCluster::setsheetid
void setsheetid(int sheetid)
Definition RecCgemCluster.h:38
RecCgemCluster::setlayerid
void setlayerid(int layerid)
Definition RecCgemCluster.h:37
RecCgemCluster::setRecZ_mTPC
void setRecZ_mTPC(double recZ)
Definition RecCgemCluster.h:49
RecCgemCluster::setSlope_mTPC
void setSlope_mTPC(double s)
Definition RecCgemCluster.h:50
RecCgemCluster::setrecv_CC
void setrecv_CC(double recv)
Definition RecCgemCluster.h:45
RecCgemCluster::setInter_mTPC
void setInter_mTPC(double q)
Definition RecCgemCluster.h:51
RecCgemCluster::getenergydeposit
double getenergydeposit(void) const
Definition RecCgemCluster.h:59
RecCgemCluster::setrecphi_CC
void setrecphi_CC(double recphi)
Definition RecCgemCluster.h:44
RecCgemCluster::setclusterid
void setclusterid(int clusterid)
Definition RecCgemCluster.h:35
RecCgemCluster::setenergydeposit
void setenergydeposit(double energydeposit)
Definition RecCgemCluster.h:40
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition RecCgemCluster.h:54
RecCgemCluster::setrecv
void setrecv(double recv)
Definition RecCgemCluster.h:42
RecCgemCluster::setRecZ
void setRecZ(double recZ)
Definition RecCgemCluster.h:43
RecCgemCluster::getlayerid
int getlayerid(void) const
Definition RecCgemCluster.h:56
RecCgemCluster::getclusterflagb
int getclusterflagb(void) const
Definition RecCgemCluster.h:73
RecCgemCluster::setrecphi_mTPC
void setrecphi_mTPC(double recphi)
Definition RecCgemCluster.h:47
RecCgemCluster::setRecZ_CC
void setRecZ_CC(double recZ)
Definition RecCgemCluster.h:46
RecCgemCluster::getrecphi
double getrecphi(void) const
Definition RecCgemCluster.h:60
RecCgemCluster::getrecv
double getrecv(void) const
Definition RecCgemCluster.h:61
RecCgemCluster::getsheetid
int getsheetid(void) const
Definition RecCgemCluster.h:57
RecCgemCluster::setflag
void setflag(int flag)
Definition RecCgemCluster.h:39
RecCgemCluster::setclusterflag
void setclusterflag(int begin, int end)
Definition RecCgemCluster.h:72
RecCgemCluster::setrecv_mTPC
void setrecv_mTPC(double recv)
Definition RecCgemCluster.h:48
RecCgemCluster::setrecphi
void setrecphi(double recphi)
Definition RecCgemCluster.h:41
RecCgemCluster::getclusterflage
int getclusterflage(void) const
Definition RecCgemCluster.h:74
Event
Definition Event.h:21