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