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