CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemLineFit.cxx
Go to the documentation of this file.
1#include "CgemLineFit/CgemLineFit.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/AlgFactory.h"
4#include "GaudiKernel/ISvcLocator.h"
5#include "GaudiKernel/SmartDataPtr.h"
6#include "GaudiKernel/IDataProviderSvc.h"
7#include "GaudiKernel/IPartPropSvc.h"
8#include "GaudiKernel/PropertyMgr.h"
9#include "GaudiKernel/IService.h"
10#include "GaudiKernel/NTuple.h"
11#include "GaudiKernel/INTupleSvc.h"
12#include "EventModel/EventHeader.h"
13#include "MdcRawEvent/MdcDigi.h"
14#include "Identifier/Identifier.h"
15#include "Identifier/MdcID.h"
16#include <iostream>
17#include <math.h>
18#include "TStopwatch.h"
19
20
21#include "EvTimeEvent/RecEsTime.h"
22#include "MdcGeom/EntranceAngle.h"
23#include "RawEvent/RawDataUtil.h"
24#include "MdcData/MdcHit.h"
25
26#include "McTruth/DecayMode.h"
27#include "McTruth/McParticle.h"
28#include "TrackUtil/Helix.h"
29#include "MdcRecoUtil/Pdt.h"
30
31#include "TrkBase/TrkFit.h"
32#include "TrkBase/TrkHitList.h"
33#include "TrkBase/TrkExchangePar.h"
34#include "TrkFitter/TrkHelixMaker.h"
35#include "TrkFitter/TrkCircleMaker.h"
36#include "TrkFitter/TrkLineMaker.h"
37#include "TrkFitter/TrkHelixFitter.h"
38#include "MdcxReco/Mdcxprobab.h"
39#include "MdcData/MdcRecoHitOnTrack.h"
40#include "MdcPrintSvc/IMdcPrintSvc.h"
41
42#include <math.h>
43
44#include "EvTimeEvent/RecEsTime.h"
45#include "McTruth/MdcMcHit.h"
46#include "TrkBase/TrkFit.h"
47#include "TrkBase/TrkHitList.h"
48#include "TrkBase/TrkExchangePar.h"
49#include "TrkFitter/TrkHelixMaker.h"
50#include "TrkFitter/TrkCircleMaker.h"
51#include "TrkFitter/TrkLineMaker.h"
52#include "TrkFitter/TrkHelixFitter.h"
53#include "MdcPrintSvc/IMdcPrintSvc.h"
54#include "MdcGeom/EntranceAngle.h"
55#include "TrackUtil/Helix.h"
56#include "GaudiKernel/IPartPropSvc.h"
57#include "MdcRecoUtil/Pdt.h"
58#include "RawEvent/RawDataUtil.h"
59#include "MdcData/MdcHit.h"
60#include "MdcTrkRecon/MdcMap.h"
61
62#include "TTree.h"
63#include "TH2D.h"
64
65#include "BesTimerSvc/IBesTimerSvc.h"
66#include "BesTimerSvc/BesTimerSvc.h"
67#include "BesTimerSvc/BesTimer.h"
68#include "TRandom.h"
69#include "TArrayI.h"
70#include "TGraph.h"
71#include "TF1.h"
72#include "TMinuit.h"
73#include "TMath.h"
74
75#include "CLHEP/Vector/ThreeVector.h"
76#include "CLHEP/Vector/LorentzVector.h"//
77#include "CLHEP/Vector/TwoVector.h"
78using CLHEP::Hep3Vector;
79using CLHEP::Hep2Vector;
80using CLHEP::HepLorentzVector;
81using namespace std;
82using namespace Event;
83
84
85vector<double>Vec_layer,Vec_phi,Vec_Z,Vec_layerid,Vec_Virlayerid,Vec_v,Vec_flag,Vec_Q,Vec_sheetid;//
86vector<double>Vec_m_layer,Vec_m_phi,Vec_m_Z,Vec_m_layerid,Vec_m_v,Vec_m_flag,Vec_m_Q,Vec_m_sheetid;//
87
88CgemMidDriftPlane* Mp;
89CgemGeoAlign* Al;
90CgemGeoLayer* GeoLayer0;
91CgemGeoLayer* GeoLayer1;
92CgemGeoLayer* GeoLayer2;
93StraightLine* l_outer;
94
95CgemGeoReadoutPlane* pl_00;
96CgemGeoReadoutPlane* pl_10;
97CgemGeoReadoutPlane* pl_11;
98CgemGeoReadoutPlane* pl_20;
99CgemGeoReadoutPlane* pl_21;
100
101bool Align_FLAG(false);
102double sigma2(0);// resoulution of hardware
103// flag for each sheet ?
104int f21(1),f11(2),f00(4),f01(8),f10(16),f20(32),fa(63),sheet_flag(0);
105
106double length,dr_set,phi0_set,dz_set,tanl_set;
107double R_layer[3];
108int NC;
109int _clst_0(0),_clst_1(0),_clst_2(0);
110
111CgemLineFit::CgemLineFit(const std::string& name, ISvcLocator* pSvcLocator) :
112 Algorithm(name, pSvcLocator)
113{
114 declareProperty("sigma", sigma = 0.13);// unit:mm
115 declareProperty("Run10_Flag", Run10_Flag = false);// true:4 clusters, false: 3 clusters
116 declareProperty("Align_Flag", Align_Flag = false);// Alignment
117 declareProperty("Switch_CCmTPC", Switch_CCmTPC = 0);// 0:Center charge 1:mTPC
118 declareProperty("Method", Method = 0); // 0: ToyMC 1: max_charge 2:closest to intersection 3: Loop_all 4: Loop_MaxQ
119 declareProperty("MAX_COMB", MAX_COMB = 50); // The maximum number of combinations for method 3
120 declareProperty("MAX_Clusters", MAX_Clusters = 500); // The maximum number of combinations for method 3
121 declareProperty("MaxClu_Sheet", Nmax = 3); //Max number of clusters on each sheet for method 4
122 declareProperty("Chisq_cut", Chisq_cut = 2000); //Max chisq for the selected set of clusters
123 declareProperty("TEST_N", TEST_N = 0); //set the sheet you want to study the efficiency / resolution. 1 to 4 from top to bottom. 0 use all the layers
124 declareProperty("Debug", debug = 0); //Switch of debug
125 declareProperty("MinQ_Clus2D", MinQ_Clus2D = 0); //Add charge cut on 2D cluster to reduce the noise
126 declareProperty("MC", MC = 0); //remove the hits on 3rd layer for MC
127}
128
129// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
130StatusCode CgemLineFit::initialize(){
131 MsgStream log(msgSvc(), name());
132 sigma2=sigma*sigma;
133 Align_FLAG=Align_Flag;
134 log << MSG::INFO << "in initialize()" << endreq;
135 StatusCode sc;
136 NTuplePtr cosmic(ntupleSvc(), "LineFit/cosmic");
137 if ( cosmic ) m_tuple = cosmic;
138 else {
139 m_tuple = ntupleSvc()->book ("LineFit/cosmic", CLID_ColumnWiseTuple, "N-Tuple example");
140 if ( m_tuple ) {
141
142 sc = m_tuple->addItem ("run",run);
143 sc = m_tuple->addItem ("event",event);
144 sc = m_tuple->addItem ("R0",R0);
145 sc = m_tuple->addItem ("R1",R1);
146 sc = m_tuple->addItem ("R2",R2);
147 sc = m_tuple->addItem ("clst_0",clst_0);
148 sc = m_tuple->addItem ("clst_1",clst_1);
149 sc = m_tuple->addItem ("clst_2",clst_2);
150 sc = m_tuple->addItem ("clst_01",clst_01);
151 sc = m_tuple->addItem ("clst_00",clst_00);
152 sc = m_tuple->addItem ("clst_11",clst_11);
153 sc = m_tuple->addItem ("clst_10",clst_10);
154 sc = m_tuple->addItem ("status1",status1);
155 sc = m_tuple->addItem ("status2",status2);
156 sc = m_tuple->addItem ("status3",status3);
157 sc = m_tuple->addItem ("ini_0",ini_0);
158 sc = m_tuple->addItem ("ini_1",ini_1);
159 sc = m_tuple->addItem ("ini_2",ini_2);
160 sc = m_tuple->addItem ("ini_3",ini_3);
161 sc = m_tuple->addItem ("inim_0",inim_0);
162 sc = m_tuple->addItem ("inim_1",inim_1);
163 sc = m_tuple->addItem ("inim_2",inim_2);
164 sc = m_tuple->addItem ("inim_3",inim_3);
165 sc = m_tuple->addItem ("inii_0",inii_0);
166 sc = m_tuple->addItem ("inii_1",inii_1);
167 sc = m_tuple->addItem ("inii_2",inii_2);
168 sc = m_tuple->addItem ("inii_3",inii_3);
169 sc = m_tuple->addItem ("par0",par0);
170 sc = m_tuple->addItem ("par1",par1);
171 sc = m_tuple->addItem ("par2",par2);
172 sc = m_tuple->addItem ("par3",par3);
173 sc = m_tuple->addItem ("MC_par0",MC_par0);
174 sc = m_tuple->addItem ("MC_par1",MC_par1);
175 sc = m_tuple->addItem ("MC_par2",MC_par2);
176 sc = m_tuple->addItem ("MC_par3",MC_par3);
177 sc = m_tuple->addItem ("MC_px",MC_px);
178 sc = m_tuple->addItem ("MC_py",MC_py);
179 sc = m_tuple->addItem ("MC_pz",MC_pz);
180 sc = m_tuple->addItem ("Err_par0",Err_par0);
181 sc = m_tuple->addItem ("Err_par1",Err_par1);
182 sc = m_tuple->addItem ("Err_par2",Err_par2);
183 sc = m_tuple->addItem ("Err_par3",Err_par3);
184
185 sc = m_tuple->addItem ("x_2_up",x_2_up);
186 sc = m_tuple->addItem ("v_2_up",v_2_up);
187 sc = m_tuple->addItem ("x_1_up",x_1_up);
188 sc = m_tuple->addItem ("v_1_up",v_1_up);
189 sc = m_tuple->addItem ("x_0_up",x_0_up);
190 sc = m_tuple->addItem ("v_0_up",v_0_up);
191 sc = m_tuple->addItem ("x_2_down",x_2_down);
192 sc = m_tuple->addItem ("v_2_down",v_2_down);
193 sc = m_tuple->addItem ("z_1_up",z_1_up);
194 sc = m_tuple->addItem ("z_2_up",z_2_up);
195 sc = m_tuple->addItem ("z_0_up",z_0_up);
196 sc = m_tuple->addItem ("z_1_down",z_1_down);
197 sc = m_tuple->addItem ("z_2_down",z_2_down);
198 sc = m_tuple->addItem ("z_0_down",z_0_down);
199 sc = m_tuple->addItem ("x_1_down",x_1_down);
200 sc = m_tuple->addItem ("v_1_down",v_1_down);
201 sc = m_tuple->addItem ("x_0_down",x_0_down);
202 sc = m_tuple->addItem ("v_0_down",v_0_down);
203
204
205 sc = m_tuple->addItem ("x_2_up_m",x_2_up_m);
206 sc = m_tuple->addItem ("x_1_up_m",x_1_up_m);
207 sc = m_tuple->addItem ("x_0_up_m",x_0_up_m);
208 sc = m_tuple->addItem ("x_0_down_m",x_0_down_m);
209 sc = m_tuple->addItem ("x_1_down_m",x_1_down_m);
210 sc = m_tuple->addItem ("x_2_down_m",x_2_down_m);
211 sc = m_tuple->addItem ("v_2_up_m",v_2_up_m);
212 sc = m_tuple->addItem ("v_1_up_m",v_1_up_m);
213 sc = m_tuple->addItem ("v_0_up_m",v_0_up_m);
214 sc = m_tuple->addItem ("v_0_down_m",v_0_down_m);
215 sc = m_tuple->addItem ("v_1_down_m",v_1_down_m);
216 sc = m_tuple->addItem ("v_2_down_m",v_2_down_m);
217 sc = m_tuple->addItem ("z_2_up_m",z_2_up_m);
218 sc = m_tuple->addItem ("z_1_up_m",z_1_up_m);
219 sc = m_tuple->addItem ("z_0_up_m",z_0_up_m);
220 sc = m_tuple->addItem ("z_0_down_m",z_0_down_m);
221 sc = m_tuple->addItem ("z_1_down_m",z_1_down_m);
222 sc = m_tuple->addItem ("z_2_down_m",z_2_down_m);
223
224 sc = m_tuple->addItem ("x_2_up_f",x_2_up_f);
225 sc = m_tuple->addItem ("x_1_up_f",x_1_up_f);
226 sc = m_tuple->addItem ("x_0_up_f",x_0_up_f);
227 sc = m_tuple->addItem ("x_0_down_f",x_0_down_f);
228 sc = m_tuple->addItem ("x_1_down_f",x_1_down_f);
229 sc = m_tuple->addItem ("x_2_down_f",x_2_down_f);
230 sc = m_tuple->addItem ("v_2_up_f",v_2_up_f);
231 sc = m_tuple->addItem ("v_1_up_f",v_1_up_f);
232 sc = m_tuple->addItem ("v_0_up_f",v_0_up_f);
233 sc = m_tuple->addItem ("v_0_down_f",v_0_down_f);
234 sc = m_tuple->addItem ("v_1_down_f",v_1_down_f);
235 sc = m_tuple->addItem ("v_2_down_f",v_2_down_f);
236
237 sc = m_tuple->addItem ("x_2_up_mc",x_2_up_mc);
238 sc = m_tuple->addItem ("x_1_up_mc",x_1_up_mc);
239 sc = m_tuple->addItem ("x_0_up_mc",x_0_up_mc);
240 sc = m_tuple->addItem ("x_0_down_mc",x_0_down_mc);
241 sc = m_tuple->addItem ("x_1_down_mc",x_1_down_mc);
242 sc = m_tuple->addItem ("x_2_down_mc",x_2_down_mc);
243 sc = m_tuple->addItem ("v_2_up_mc",v_2_up_mc);
244 sc = m_tuple->addItem ("v_1_up_mc",v_1_up_mc);
245 sc = m_tuple->addItem ("v_0_up_mc",v_0_up_mc);
246 sc = m_tuple->addItem ("v_0_down_mc",v_0_down_mc);
247 sc = m_tuple->addItem ("v_1_down_mc",v_1_down_mc);
248 sc = m_tuple->addItem ("v_2_down_mc",v_2_down_mc);
249
250
251 sc = m_tuple->addItem ("inter_1_x",inter_1_x);
252 sc = m_tuple->addItem ("inter_2_x",inter_2_x);
253 sc = m_tuple->addItem ("inter_3_x",inter_3_x);
254 sc = m_tuple->addItem ("inter_4_x",inter_4_x);
255
256 sc = m_tuple->addItem ("inter_1_V",inter_1_V);
257 sc = m_tuple->addItem ("inter_2_V",inter_2_V);
258 sc = m_tuple->addItem ("inter_3_V",inter_3_V);
259 sc = m_tuple->addItem ("inter_4_V",inter_4_V);
260
261 sc = m_tuple->addItem ("inter_1_flag",inter_1_flag);
262 sc = m_tuple->addItem ("inter_2_flag",inter_2_flag);
263 sc = m_tuple->addItem ("inter_3_flag",inter_3_flag);
264 sc = m_tuple->addItem ("inter_4_flag",inter_4_flag);
265
266 sc = m_tuple->addItem ("inter_1_z",inter_1_z);
267 sc = m_tuple->addItem ("inter_2_z",inter_2_z);
268 sc = m_tuple->addItem ("inter_3_z",inter_3_z);
269 sc = m_tuple->addItem ("inter_4_z",inter_4_z);
270
271
272 sc = m_tuple->addItem ("chisq_1",chisq_1);
273 sc = m_tuple->addItem ("chisq_2",chisq_2);
274 sc = m_tuple->addItem ("chisq_3",chisq_3);
275 sc = m_tuple->addItem ("pos_x",pos_x);
276 sc = m_tuple->addItem ("pos_y",pos_y);
277 sc = m_tuple->addItem ("pos_z",pos_z);
278
279 sc = m_tuple->addItem ("hit_x",hit_x);
280 sc = m_tuple->addItem ("hit_y",hit_y);
281 sc = m_tuple->addItem ("hit_z",hit_z);
282
283 sc = m_tuple->addItem ("ang_1",ang_1);
284 sc = m_tuple->addItem ("ang_1_x",ang_1_x);
285 sc = m_tuple->addItem ("ang_1_v",ang_1_v);
286
287 sc = m_tuple->addItem ("ang_2",ang_2);
288 sc = m_tuple->addItem ("ang_2_x",ang_2_x);
289 sc = m_tuple->addItem ("ang_2_v",ang_2_v);
290
291 sc = m_tuple->addItem ("ang_3",ang_3);
292 sc = m_tuple->addItem ("ang_3_x",ang_3_x);
293 sc = m_tuple->addItem ("ang_3_v",ang_3_v);
294
295 sc = m_tuple->addItem ("ang_4",ang_4);
296 sc = m_tuple->addItem ("ang_4_x",ang_4_x);
297 sc = m_tuple->addItem ("ang_4_v",ang_4_v);
298 sc = m_tuple->addItem ("Res_phi",Res_phi);
299 sc = m_tuple->addItem ("Res_V",Res_V);
300 sc = m_tuple->addItem ("Test_phi",Test_phi);
301 sc = m_tuple->addItem ("Test_V",Test_V);
302 }
303 else {
304 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple) << endmsg;
305 return StatusCode::FAILURE;
306 }
307 }
308
309
310
311 NTuplePtr Track(ntupleSvc(), "LineFit/Track");
312 if ( Track ) m_tuple_track = Track;
313 else {
314 m_tuple_track = ntupleSvc()->book ("LineFit/Track", CLID_ColumnWiseTuple, "N-Tuple example");
315 if ( m_tuple_track ) {
316 sc = m_tuple_track->addItem ("run", tr_run);
317 sc = m_tuple_track->addItem ("event", tr_event);
318 sc = m_tuple_track->addItem ("drho", tr_drho);
319 sc = m_tuple_track->addItem ("phi0", tr_phi0);
320 sc = m_tuple_track->addItem ("dz0", tr_dz0);
321 sc = m_tuple_track->addItem ("chisq", tr_chisq);
322 sc = m_tuple_track->addItem ("Is_fitted", tr_Is_fitted);
323 sc = m_tuple_track->addItem ("tanLam", tr_tanLam);
324 sc = m_tuple_track->addItem ("ncluster", tr_ncluster, 0, 5000);
325 sc = m_tuple_track->addItem ("id_cluster", tr_ncluster, tr_id_cluster);
326 sc = m_tuple_track->addItem ("x_cluster", tr_ncluster, tr_x_cluster);
327 sc = m_tuple_track->addItem ("y_cluster", tr_ncluster, tr_y_cluster);
328 sc = m_tuple_track->addItem ("z_cluster", tr_ncluster, tr_z_cluster);
329 sc = m_tuple_track->addItem ("Q_cluster", tr_ncluster, tr_Q_cluster);
330 sc = m_tuple_track->addItem ("phi_cluster", tr_ncluster, tr_phi_cluster);
331 sc = m_tuple_track->addItem ("layer_cluster", tr_ncluster, tr_layer_cluster);
332 sc = m_tuple_track->addItem ("sheet_cluster", tr_ncluster, tr_sheet_cluster);
333 sc = m_tuple_track->addItem ("Is_selected_cluster", tr_ncluster, tr_Is_selected_cluster);
334
335 }
336 else {
337 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple_track) << endmsg;
338 return StatusCode::FAILURE;
339 }
340 }
341
342 IRawDataProviderSvc* irawDataProviderSvc;
343 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
344 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
345 if ( sc.isFailure() ){
346 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
347 return StatusCode::FAILURE;
348 }
349
350 ICgemGeomSvc* icgemGeomSvc;
351 sc = service ("CgemGeomSvc", icgemGeomSvc);
352 m_cgemGeomSvc = dynamic_cast<CgemGeomSvc*> (icgemGeomSvc);
353
354 if ( sc.isFailure() ){
355 log << MSG::FATAL << "Could not load CgemGeomSvc!" << endreq;
356 return StatusCode::FAILURE;
357 }
358
359 R_layer[0]=(m_cgemGeomSvc->getCgemLayer(0)->getMiddleROfGapD());
360 R_layer[1]=(m_cgemGeomSvc->getCgemLayer(1)->getMiddleROfGapD());
361 R_layer[2]=(m_cgemGeomSvc->getCgemLayer(2)->getMiddleROfGapD());
362 length= m_cgemGeomSvc->getLengthOfCgem();
363 Mp=m_cgemGeomSvc->getMidDriftPtr() ;
364 Al=m_cgemGeomSvc->getAlignPtr() ;
365
366 GeoLayer0=m_cgemGeomSvc->getCgemLayer(0);
367 GeoLayer1=m_cgemGeomSvc->getCgemLayer(1);
368 GeoLayer2=m_cgemGeomSvc->getCgemLayer(2);
369
370 pl_00 =m_cgemGeomSvc->getReadoutPlane( 0, 0) ;
371 pl_10=m_cgemGeomSvc->getReadoutPlane( 1, 0) ;
372 pl_11 =m_cgemGeomSvc->getReadoutPlane( 1, 1) ;
373 pl_20=m_cgemGeomSvc->getReadoutPlane( 2, 0) ;
374 pl_21 =m_cgemGeomSvc->getReadoutPlane( 2, 1) ;
375
376
377 cout << "CgemLineFit alignment: is it on? " << Align_Flag << endl;
378 if(Align_Flag) {
379 for(int ilay=0; ilay<6; ilay++) {
380 cout << "LAYER " << ilay+1 << endl;
381 cout << "shift dx " << Al->getDx(ilay) << " dy " << Al->getDy(ilay) << " dz " << Al->getDz(ilay) << endl;
382 cout << "rotation Rx " << Al->getRx(ilay) << " Ry " << Al->getRy(ilay) << " Rz " << Al->getRz(ilay) << endl;
383 cout << "midplane radius " << Mp->getR(int(ilay/2)) << endl;
384 }
385 }
386
387 ICgemCalibFunSvc* icgemCalibSvc;
388 sc = service ("CgemCalibFunSvc", icgemCalibSvc);
389 m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc*>(icgemCalibSvc);
390 if ( sc.isFailure() ){
391 log << MSG::FATAL << "Could not load CgemCalibFunSvc!" << endreq;
392 return StatusCode::FAILURE;
393 }
394
395 return StatusCode::SUCCESS;
396}
397
398// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
399StatusCode CgemLineFit::execute() {
400 MsgStream log(msgSvc(), name());
401 log << MSG::INFO << "in execute()" << endreq;
402
403 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
404 if (!eventHeader) {
405 log << MSG::FATAL << "Could not find Event Header" << endreq;
406 return StatusCode::FAILURE;
407 }
408
409 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
410 if (!recCgemClusterCol){
411 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
412 return StatusCode::FAILURE;
413 }
414
415 m_event=eventHeader->eventNumber();
416 m_run=eventHeader->runNumber();
417 event=m_event;
418 run=m_run;
419 tr_run=run;
420 tr_event=event;
421 if((m_event % 100)==0)cout<<" begin evt : "<<event <<endl;
422 //////////////////Collecting Clusters///////////
423
424 // release the memory
425 Vec_layer.clear();
426 Vec_phi.clear();
427 Vec_Z.clear();
428 Vec_v.clear();
429 Vec_layerid.clear();
430 Vec_Virlayerid.clear();
431 Vec_flag.clear();
432 Vec_Q.clear();
433 Vec_sheetid.clear();
434 Vec_m_layer.clear();
435 Vec_m_phi.clear();
436 Vec_m_Z.clear();
437 Vec_m_v.clear();
438 Vec_m_layerid.clear();
439 Vec_m_flag.clear();
440 Vec_m_Q.clear();
441 Vec_m_sheetid.clear();
442 vecclusters.clear();
443 SmartRefVector<RecCgemCluster>().swap(vecclusters);
444
445 _clst_0=0;_clst_1=0;_clst_2=0;
446
447 NXComb = -1;
448 NVComb = -1;
449
450 TStopwatch gtimer;
451 gtimer.Start();
452
453 if(debug) cout<<"start Straight line fit"<<endl;
454 if(Method==0){
455 if(!ToyMC()){
456 tr_Is_fitted = 0;
457 tr_chisq = -1;
458 m_tuple_track->write();
459 return StatusCode::SUCCESS;}}
460 else if(Method==1){
461 if(!Data_Max()){
462 tr_Is_fitted = 0;
463 tr_chisq = -1;
464 m_tuple_track->write();
465 return StatusCode::SUCCESS;}}
466 else if(Method==2){
467 if(!Data_Closest()){
468 tr_Is_fitted = 0;
469 tr_chisq = -1;
470 m_tuple_track->write();
471 return StatusCode::SUCCESS;}}
472 else if(Method==3){
473 if(!Loop_All()){
474 tr_Is_fitted = 0;
475 tr_chisq = -1;
476 m_tuple_track->write();
477 return StatusCode::SUCCESS;}}
478 else if(Method==4){
479 if(!Loop_MaxQ()){
480 tr_Is_fitted = 0;
481 tr_chisq = -1;
482 m_tuple_track->write();
483 return StatusCode::SUCCESS;}}
484
485 if(Vec_layerid.size()<4&&Run10_Flag)
486 return StatusCode::SUCCESS;
487
488 /////Set initial parameter according to the clusters on the outtermost layer and perform the final fit:
489 /////1, fit phi0, drho; 2, fit tanl z0; 3, fit 4 par. at the same time/////
490
491 TMinuit* fit=Fit(l_outer,fa,0);
492 double trkpar[4]={-999,-999,-999,-999};
493 double trkparErr[4]={-999,-999,-999,-999};
494 Int_t ierflg ,npari,nparx,istat3;
495 Double_t fmin,edm,errdef;
496 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat3);
497 for(int i=0; i<4; i++){
498 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
499 double emat[16];
500 fit->mnemat(emat,4);
501 chisq_3=fit->fAmin;
502 registerTrack(m_trackList_tds,m_hitList_tds);
503 Store_Trk(fit,0);
504 delete fit;
505 double chisq_last = chisq_3;
506
507 gtimer.Stop();
508 vector<double> _trkpar=Trans(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
509
510 par0=_trkpar[0];
511 par1=_trkpar[1];
512 par2=_trkpar[2];
513 par3=_trkpar[3];
514 Err_par0=trkparErr[0];
515 Err_par1=trkparErr[1];
516 Err_par2=trkparErr[2];
517 Err_par3=trkparErr[3];
518
519 status3=istat3;
520
521 clst_0=_clst_0;
522 clst_1=_clst_1;
523 clst_2=_clst_2;
524
525 R0=R_layer[0];
526 R1=R_layer[1];
527 R2=R_layer[2];
528
529 // the initial track parameters based on the clusters on the outtermost layer
530 ini_0=l_outer->dr();
531 ini_1=l_outer->phi0();
532 ini_2=l_outer->dz();
533 ini_3=l_outer->tanl();
534
535 if(debug) cout<<"Get_otherIniPar"<<endl;
536 // Get the initial track parameters based on the clusters on middle and innermost layers
537 if(TEST_N==0) Get_OtherIniPar(_clst_0,_clst_1,_clst_2);
538
539 // Store the x,v,z information of clusters on the avialable layers,_clst_0: N clusters on middle layer, _clst_1: N clusters on innermost layer
540 // _clst_0=2;
541 // _clst_1=2;
542 // _clst_2=0;
543
544 if(debug) cout<<"Rec_Clusters"<<endl;
545 Rec_Clusters();
546 // Store the x,v,z information of clusters by mTPC method on the avialable layers
547 if(Switch_CCmTPC==1) Rec_Clusters_mTPC();
548
549 if(debug) cout<<"Fit_Clusters"<<endl;
550 // Store the intersections between fitted line and Cgem
551 Fit_Clusters(trkpar);
552
553 if(debug) cout<<"GetIntersection"<<endl;
554 // Store the intersections according to the extroplation of the clusters on other layers
555 GetIntersection();
556
557 if(debug) cout<<"GetRes"<<endl;
558 // Get the resolution
559 StraightLine l_fit(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
560 GetRes(&l_fit,TEST_N);
561 m_tuple->write();
562
563 tr_drho =trkpar[0];
564 tr_phi0 =trkpar[1];
565 tr_dz0 =trkpar[2];
566 tr_tanLam=trkpar[3];
567 tr_Is_fitted=1;
568 tr_chisq=chisq_last;
569 m_tuple_track->write();
570 if(debug) cout<<"Finish all"<<endl;
571
572 delete l_outer;
573 return StatusCode::SUCCESS;
574
575}
576
577StatusCode CgemLineFit::finalize() {
578 MsgStream log(msgSvc(), name());
579 log << MSG::INFO<< "in finalize()" << endreq;
580 return StatusCode::SUCCESS;
581}
582
583//========================================================================================================
584// Method0: for ToyMC study
585//========================================================================================================
586
587bool CgemLineFit::ToyMC() {
588
589 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
590 int _loop(0);
591 if(!Get_MCInf()) return false;
592 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
593 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
594 for(;iter!=recCgemClusterCol->end();iter++)
595 {
596 _loop++;
597 RecCgemCluster *cluster = (*iter);
598 int flag=cluster->getflag();
599 if(flag!=2) continue;
600
601 int layerid=cluster->getlayerid();
602 int sheetid=cluster->getsheetid();
603 double _phi=cluster->getrecphi();
604 if(MC) _phi = RealPhi(_phi);
605 double _v=cluster->getrecv();
606 double _Z=cluster->getRecZ();
607 double _Q=cluster->getenergydeposit();
608 int vir_layerid = GetVirLay(layerid, _phi);
609 if(layerid==2||((!Run10_Flag)&&layerid==0&&_phi>acos(-1)))continue;
610
611 Vec_flag.push_back(flag);
612 Vec_layerid.push_back(layerid);
613 Vec_Virlayerid.push_back(vir_layerid);
614 Vec_sheetid.push_back(sheetid);
615 Vec_layer.push_back(R_layer[layerid]);
616
617 Vec_phi.push_back(_phi);
618 Vec_v.push_back(_v );
619 Vec_Z.push_back(_Z);
620
621 Vec_Q.push_back(_Q);
622 cluster->setTrkId(0);
623 vecclusters.push_back(cluster);
624
625 }
626 int NC = vecclusters.size();
627 if(debug) cout<<"NC is "<<NC<<endl;
628
629 if(Vec_layer.size()<3)
630 {
631 return false;
632 }
633
634 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
635 {
636 if(Vec_layerid[i_cl]==2) _clst_2++;
637 if(Vec_layerid[i_cl]==1) _clst_1++;
638 if(Vec_layerid[i_cl]==0) _clst_0++;
639 }
640
641 if(debug) cout<<"clst1, clst0 are "<<_clst_1<<", "<<_clst_0<<endl;
642 OrderClusters();
643 if(_clst_1==2)
644 l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2); // should be this
645 else return false;
646
647 if(_clst_0>2)Filter(0,l_outer);
648
649 return true;
650}
651
652//========================================================================================================
653// Method1: select the cluster with the largest charge on each hemisphere
654//========================================================================================================
655
656bool CgemLineFit::Data_Max() {
657 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
658 int _loop(0);
659 double maxQ_11(0),maxQ_10(0),maxQ_00(0),maxQ_01(0);
660 int max_11(-1),max_10(-1),max_00(-1),max_01(-1);
661 double maxv_00(0),maxv_01;
662
663 double phi_11(-1),phi_10(-1),phi_00(-1),phi_01(-1);
664 double z_11(0),z_10(0),z_00(0),z_01(0);
665 double cid_11(-1),cid_10(-1),cid_00(-1),cid_01(-1);
666 double Xid_11(-1),Xid_10(-1),Xid_00(-1),Xid_01(-1);
667 double Vid_11(-1),Vid_10(-1),Vid_00(-1),Vid_01(-1);
668
669 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
670 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
671
672 int ic(0);
673 int ic_tot(0);
674 for(;iter!=recCgemClusterCol->end();iter++)
675 {
676 _loop++;
677 RecCgemCluster *cluster = (*iter);
678 int flag=cluster->getflag();
679 int layerid=cluster->getlayerid();
680 double _Q=cluster->getenergydeposit();
681 double _phi=cluster->getrecphi();
682 double _v=cluster->getrecv();
683 int sheetid=cluster->getsheetid();
684
685 int clusterid=cluster->getclusterid();
686 int clusterflagb=cluster->getclusterflagb();
687 int clusterflage=cluster->getclusterflage();
688 double x=R_layer[layerid]*cos(_phi);
689 double y=R_layer[layerid]*sin(_phi);
690 double z=cluster->getRecZ();
691 if(_loop>MAX_Clusters) break;
692
693 if(flag!=2)continue;
694 if(_Q<MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
695
696 tr_id_cluster[ic]=clusterid;
697 tr_x_cluster[ic]=x;
698 tr_y_cluster[ic]=y;
699 tr_z_cluster[ic]=z;
700 tr_Q_cluster[ic]=_Q;
701 tr_phi_cluster[ic]=_phi;
702 tr_layer_cluster[ic]=layerid;
703 tr_sheet_cluster[ic]=sheetid;
704
705 ic_tot++;
706 ic++;
707 if(layerid==1&&sheetid==1&&maxQ_11<_Q){
708 maxQ_11=_Q;max_11=_loop;phi_11=_phi;z_11=z;cid_11=clusterid;Xid_11=clusterflagb;Vid_11=clusterflage;}
709 if(layerid==1&&sheetid==0&&maxQ_10<_Q){
710 maxQ_10=_Q;max_10=_loop;phi_10=_phi;z_10=z;cid_10=clusterid;Xid_10=clusterflagb;Vid_10=clusterflage;}
711 if(layerid==0&&sheetid==0&&_phi>0&&maxQ_01<_Q&&maxv_00!=_v){
712 maxQ_01=_Q;max_01=_loop;phi_01=_phi;z_01=z;cid_01=clusterid;Xid_01=clusterflagb;Vid_01=clusterflage;}
713 if(layerid==0&&sheetid==0&&_phi<0&&maxQ_00<_Q&&maxv_00!=_v){
714 maxQ_00=_Q;max_00=_loop;phi_00=_phi;z_00=z;cid_00=clusterid;Xid_00=clusterflagb;Vid_00=clusterflage;}
715
716 if(layerid==1&&sheetid==1)cl_11++;
717 if(layerid==0&&sheetid==0&&_phi>0)cl_01++;
718 if(layerid==0&&sheetid==0&&_phi<0)cl_00++;
719 if(layerid==1&&sheetid==0)cl_10++;
720 }
721
722
723 tr_ncluster = ic;
724 for(int i=0; i<ic; i++)
725 {
726 tr_Is_selected_cluster[i]=0;
727 if(i==(max_00-1)&&maxQ_00>0) tr_Is_selected_cluster[i] =1;
728 if(i==(max_01-1)&&maxQ_01>0) tr_Is_selected_cluster[i] =1;
729 if(i==(max_10-1)&&maxQ_10>0) tr_Is_selected_cluster[i] =1;
730 if(i==(max_11-1)&&maxQ_11>0) tr_Is_selected_cluster[i] =1;
731 }
732
733
734 clst_11=cl_11;
735 clst_01=cl_01;
736 clst_10=cl_10;
737 clst_00=cl_00;
738 _loop=0;
739
740 iter =recCgemClusterCol->begin();
741 for(;iter!=recCgemClusterCol->end();iter++)
742 {
743 _loop++;
744 RecCgemCluster *cluster = (*iter);
745 int flag=cluster->getflag();
746 if(flag!=2)continue;
747
748 int layerid=cluster->getlayerid();
749 int sheetid=cluster->getsheetid();
750 double _phi=cluster->getrecphi();
751 double _v=cluster->getrecv();
752 double _Z=cluster->getRecZ();
753 double _Q=cluster->getenergydeposit();
754 double t_phi=cluster->getrecphi_mTPC();
755 double t_v=cluster->getrecv_mTPC();
756 double t_Z=cluster->getRecZ_mTPC();
757
758 if(!(_loop==max_00||_loop==max_10||_loop==max_01||_loop==max_11))continue;
759
760 Vec_flag.push_back(flag);
761 Vec_layerid.push_back(layerid);
762 Vec_sheetid.push_back(sheetid);
763 Vec_layer.push_back(R_layer[layerid]);
764
765 if(Switch_CCmTPC==1)
766 {
767 Vec_m_phi.push_back(t_phi);
768 Vec_m_v.push_back(t_v);
769 Vec_m_Z.push_back(t_Z);
770 }
771 else if(Switch_CCmTPC==0)
772 {
773 Vec_phi.push_back(_phi);
774 Vec_v.push_back(_v );
775 Vec_Z.push_back(_Z);
776
777 Vec_m_phi.push_back(t_phi);
778 Vec_m_v.push_back(t_v);
779 Vec_m_Z.push_back(t_Z);
780
781 }
782
783 Vec_Q.push_back(_Q);
784 cluster->setTrkId(0);
785 vecclusters.push_back(cluster);
786
787 }
788 NC = vecclusters.size();
789
790 if(Vec_layer.size()!=4) return false;
791 NXComb = 1;
792 NVComb = 1;
793
794 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
795 {
796 if(Vec_layerid[i_cl]==2)_clst_2++;
797 if(Vec_layerid[i_cl]==1)_clst_1++;
798 if(Vec_layerid[i_cl]==0)_clst_0++;
799 }
800
801 OrderClusters();
802
803 if(_clst_1==2)
804 l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
805 else return false;
806
807 return true;
808}
809
810//========================================================================================================
811// Method2: select the cluster closest to the intersection. probably not be finished?
812//========================================================================================================
813
814bool CgemLineFit::Data_Closest() {
815 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
816 int _loop(0);
817 double maxQ_11(0),maxQ_10(0),maxQ_00(0);
818 int max_11(-1),max_10(-1),max_00(-1);
819 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
820 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
821
822 for(;iter!=recCgemClusterCol->end();iter++)
823 {
824 _loop++;
825 RecCgemCluster *cluster = (*iter);
826 int flag=cluster->getflag();
827 int layerid=cluster->getlayerid();
828 double _Q=cluster->getenergydeposit();
829 int sheetid=cluster->getsheetid();
830 double _phi=cluster->getrecphi();
831 if(flag!=2)continue;
832 if(layerid==1&&sheetid==1&&maxQ_11<_Q){
833 maxQ_11=_Q;max_11=_loop; }
834 if(layerid==1&&sheetid==0&&maxQ_10<_Q){
835 maxQ_10=_Q;max_10=_loop; }
836 if(layerid==1&&sheetid==1)cl_11++;
837 if(layerid==0&&sheetid==0&&_phi>0)cl_01++;
838 if(layerid==0&&sheetid==0&&_phi<0)cl_00++;
839 if(layerid==1&&sheetid==0)cl_10++;
840 }
841 clst_11=cl_11;
842 clst_01=cl_01;
843 clst_10=cl_10;
844 clst_00=cl_00;
845 _loop=0;
846 iter =recCgemClusterCol->begin();
847 for(;iter!=recCgemClusterCol->end();iter++)
848 {
849 _loop++;
850 RecCgemCluster *cluster = (*iter);
851 int flag=cluster->getflag();
852 if(flag!=2)continue;
853
854 int layerid=cluster->getlayerid();
855 int sheetid=cluster->getsheetid();
856 double _phi=cluster->getrecphi();
857 double _v=cluster->getrecv();
858 double _Z=cluster->getRecZ();
859 double _Q=cluster->getenergydeposit();
860 double t_phi=cluster->getrecphi_mTPC();
861 double t_v=cluster->getrecv_mTPC();
862 double t_Z=cluster->getRecZ_mTPC();
863
864 if(_loop!=max_11&&_loop!=max_10&&layerid!=0)continue;
865 Vec_flag.push_back(flag);
866 Vec_layerid.push_back(layerid);
867 Vec_sheetid.push_back(sheetid);
868 Vec_layer.push_back(R_layer[layerid]);
869 if(Switch_CCmTPC==1)
870 {
871 Vec_phi.push_back(t_phi);
872 Vec_v.push_back(t_v);
873 Vec_Z.push_back(t_Z);
874 Vec_m_phi.push_back(_phi);
875 Vec_m_v.push_back(_v);
876 Vec_m_Z.push_back(_Z);
877 }
878 else if(Switch_CCmTPC==0)
879 {
880 Vec_phi.push_back(_phi);
881 Vec_v.push_back(_v );
882 Vec_Z.push_back(_Z);
883
884 Vec_m_phi.push_back(t_phi);
885 Vec_m_v.push_back(t_v);
886 Vec_m_Z.push_back(t_Z);
887
888 }
889
890 Vec_Q.push_back(_Q);
891 cluster->setTrkId(0);
892 vecclusters.push_back(cluster);
893
894 }
895 NC = vecclusters.size();
896
897 if(Vec_layer.size()<3)
898 {
899 return false;
900 }
901 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
902 {
903 if(Vec_layerid[i_cl]==2)_clst_2++;
904 if(Vec_layerid[i_cl]==1)_clst_1++;
905 if(Vec_layerid[i_cl]==0)_clst_0++;
906 }
907
908 OrderClusters();
909
910 if(_clst_1==2)
911 l_outer=IniPar(Vec_phi[0],Vec_Z[0],1,Vec_phi[1],Vec_Z[1],1);
912 else return false;
913
914 if(!Layer_cross(0,l_outer))return false;
915 if(_clst_0>1||(Run10_Flag&&_clst_0>2))Filter(0,l_outer);
916
917 return true;
918}
919
920
921//========================================================================================================
922// Method3: loop all the clusters then select the combination giving the smallest chisq
923//========================================================================================================
924
925bool CgemLineFit::Loop_All() {
926 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
927 int _loop(0);
928 double maxQ_11(0),maxQ_10(0),maxQ_00(0);
929 int max_11(-1),max_10(-1),max_00(-1);
930 double C_00(0),C_10(0),C_11(0),C_01(0);
931 vector<double>Vec_00_layer,Vec_00_phi,Vec_00_Z,Vec_00_layerid,Vec_00_v,Vec_00_flag,Vec_00_Q,Vec_00_sheetid;//
932 vector<double>Vec_01_layer,Vec_01_phi,Vec_01_Z,Vec_01_layerid,Vec_01_v,Vec_01_flag,Vec_01_Q,Vec_01_sheetid;//
933 vector<double>Vec_10_layer,Vec_10_phi,Vec_10_Z,Vec_10_layerid,Vec_10_v,Vec_10_flag,Vec_10_Q,Vec_10_sheetid;//
934 vector<double>Vec_11_layer,Vec_11_phi,Vec_11_Z,Vec_11_layerid,Vec_11_v,Vec_11_flag,Vec_11_Q,Vec_11_sheetid;//
935
936 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
937 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
938 vector<int>L_00,L_01,L_10,L_11;
939 int ic(0);
940 for(;iter!=recCgemClusterCol->end();iter++)
941 {
942 _loop++;
943 RecCgemCluster *cluster = (*iter);
944
945 int flag=cluster->getflag();
946 int layerid=cluster->getlayerid();
947 double _Q=cluster->getenergydeposit();
948 double _phi=cluster->getrecphi();
949 double _v=cluster->getrecv();
950 int sheetid=cluster->getsheetid();
951 int clusterid=cluster->getclusterid();
952 int clusterflagb=cluster->getclusterflagb();
953 int clusterflage=cluster->getclusterflage();
954 double x=R_layer[layerid]*cos(_phi);
955 double y=R_layer[layerid]*sin(_phi);
956 double z=cluster->getRecZ();
957 if(_loop>MAX_Clusters) break;
958 ic++;
959 if(flag!=0)continue;
960
961 if(layerid==1&&sheetid==1&&TEST_N!=1){
962 L_11.push_back(_loop); }
963 if(layerid==1&&sheetid==0&&TEST_N!=4){
964 L_10.push_back(_loop); }
965 if(layerid==0&&sheetid==0&&_phi>0&&TEST_N!=2){
966 L_01.push_back(_loop); }
967 if(layerid==0&&sheetid==0&&_phi<0&&TEST_N!=3){
968 L_00.push_back(_loop); }
969 }
970
971 if(TEST_N==1) L_11.push_back(-1);
972 if(TEST_N==2) L_01.push_back(-1);
973 if(TEST_N==3) L_00.push_back(-1);
974 if(TEST_N==4) L_10.push_back(-1);
975
976 double Min_chi(1E10);
977
978 if((L_00.size()*L_01.size()*L_10.size()*L_11.size())>MAX_COMB)
979 {
980 return false;
981 }
982
983 NXComb = L_00.size()*L_01.size()*L_10.size()*L_11.size();
984
985 for(int i_11=0;i_11<L_11.size();i_11++){
986 for(int i_10=0;i_10<L_10.size();i_10++){
987 for(int i_00=0;i_00<L_00.size();i_00++){
988 for(int i_01=0;i_01<L_01.size();i_01++){
989
990 Vec_layer.clear();
991 Vec_phi.clear();
992 Vec_Z.clear();
993 Vec_v.clear();
994 Vec_layerid.clear();
995 Vec_Virlayerid.clear();
996 Vec_flag.clear();
997 Vec_Q.clear();
998 Vec_sheetid.clear();
999 Vec_m_layer.clear();
1000 Vec_m_phi.clear();
1001 Vec_m_Z.clear();
1002 Vec_m_v.clear();
1003 Vec_m_layerid.clear();
1004 Vec_m_flag.clear();
1005 Vec_m_Q.clear();
1006 Vec_m_sheetid.clear();
1007 _loop=0;
1008 iter =recCgemClusterCol->begin();
1009 for(;iter!=recCgemClusterCol->end();iter++)
1010 {
1011 _loop++;
1012 RecCgemCluster *cluster = (*iter);
1013 int flag=cluster->getflag();
1014 if(!(_loop==L_11[i_11]||_loop==L_10[i_10]||_loop==L_00[i_00]||_loop==L_01[i_01]))continue;
1015 int layerid=cluster->getlayerid();
1016 int sheetid=cluster->getsheetid();
1017 double _phi=cluster->getrecphi();
1018 double _v=cluster->getrecv();
1019 double _Z=cluster->getRecZ();
1020 double _Q=cluster->getenergydeposit();
1021 double t_phi=cluster->getrecphi_mTPC();
1022 double t_v=cluster->getrecv_mTPC();
1023 double t_Z=cluster->getRecZ_mTPC();
1024
1025 int vir_layerid = GetVirLay(layerid, _phi);
1026
1027 Vec_flag.push_back(flag);
1028 Vec_layerid.push_back(layerid);
1029 Vec_Virlayerid.push_back(vir_layerid);
1030 Vec_sheetid.push_back(sheetid);
1031 Vec_layer.push_back(R_layer[layerid]);
1032 if(Switch_CCmTPC==1)
1033 {
1034 Vec_phi.push_back(t_phi);
1035 Vec_v.push_back(t_v);
1036 Vec_Z.push_back(t_Z);
1037 Vec_m_phi.push_back(_phi);
1038 Vec_m_v.push_back(_v);
1039 Vec_m_Z.push_back(_Z);
1040 }
1041 else if(Switch_CCmTPC==0)
1042 {
1043 Vec_phi.push_back(_phi);
1044 Vec_v.push_back(_v );
1045 Vec_Z.push_back(_Z);
1046
1047 Vec_m_phi.push_back(t_phi);
1048 Vec_m_v.push_back(t_v);
1049 Vec_m_Z.push_back(t_Z);
1050
1051 }
1052 }// one combo
1053 if(Vec_layerid.size()!=4&&TEST_N==0)continue;
1054 else if(Vec_layerid.size()!=3&&TEST_N>0)continue;
1055
1056 OrderClusters();
1057
1058 StraightLine* l_outer_tmp;
1059
1060 if(TEST_N==1||TEST_N==4) l_outer_tmp=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1061 else if(TEST_N==2||TEST_N==3) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1062 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1063
1064 TMinuit* _fit=Fit(l_outer_tmp,fa,1);
1065 double _chi=_fit->fAmin;
1066 if(_chi<Min_chi)
1067 {
1068 Min_chi=_chi;
1069 if(TEST_N==0) C_00=Vec_phi[3];
1070 else C_00=-99;
1071 C_01=Vec_phi[2];
1072 C_10=Vec_phi[1];
1073 C_11=Vec_phi[0];
1074
1075 }
1076 delete _fit;
1077 delete l_outer_tmp;
1078 }
1079 }
1080 }
1081 }
1082
1083
1084 iter =recCgemClusterCol->begin();
1085 L_00.clear();L_01.clear();L_10.clear();L_11.clear();
1086 _loop=0;
1087 int ic_tot(0);
1088 for(;iter!=recCgemClusterCol->end();iter++)
1089 {
1090 _loop++;
1091 RecCgemCluster *cluster = (*iter);
1092 int flag=cluster->getflag();
1093 if(flag!=2)continue;
1094 ic_tot++;
1095 int layerid=cluster->getlayerid();
1096 double _Q=cluster->getenergydeposit();
1097 double _phi=cluster->getrecphi();
1098 int sheetid=cluster->getsheetid();
1099 if(_Q<MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
1100 if(C_00==_phi||C_01==_phi||C_11==_phi||C_10==_phi){
1101 if(layerid==1&&sheetid==1){
1102 L_11.push_back(_loop); }
1103 if(layerid==1&&sheetid==0){
1104 L_10.push_back(_loop); }
1105 if(layerid==0&&sheetid==0&&_phi>0){
1106 L_01.push_back(_loop); }
1107 if(layerid==0&&sheetid==0&&_phi<0){
1108 L_00.push_back(_loop); }
1109 }
1110
1111 }
1112 Min_chi=1E10;
1113
1114 if(TEST_N==1)L_11.push_back(-1);
1115 if(TEST_N==2)L_01.push_back(-1);
1116 if(TEST_N==3)L_00.push_back(-1);
1117 if(TEST_N==4)L_10.push_back(-1);
1118
1119 clst_11=L_11.size();
1120 clst_01=L_01.size();
1121 clst_10=L_10.size();
1122 clst_00=L_00.size();
1123 if((L_00.size()*L_01.size()*L_10.size()*L_11.size())>MAX_COMB)
1124 {cout<<"xv size much "<<endl;
1125 return false;
1126 }
1127
1128 NVComb = L_00.size()*L_01.size()*L_10.size()*L_11.size();
1129
1130 for(int i_11=0;i_11<L_11.size();i_11++){
1131 for(int i_10=0;i_10<L_10.size();i_10++){
1132 for(int i_00=0;i_00<L_00.size();i_00++){
1133 for(int i_01=0;i_01<L_01.size();i_01++){
1134
1135 Vec_layer.clear();
1136 Vec_phi.clear();
1137 Vec_Z.clear();
1138 Vec_v.clear();
1139 Vec_layerid.clear();
1140 Vec_Virlayerid.clear();
1141 Vec_flag.clear();
1142 Vec_Q.clear();
1143 Vec_sheetid.clear();
1144 Vec_m_layer.clear();
1145 Vec_m_phi.clear();
1146 Vec_m_Z.clear();
1147 Vec_m_v.clear();
1148 Vec_m_layerid.clear();
1149 Vec_m_flag.clear();
1150 Vec_m_Q.clear();
1151 Vec_m_sheetid.clear();
1152 _loop=0;
1153 iter =recCgemClusterCol->begin();
1154 for(;iter!=recCgemClusterCol->end();iter++)
1155 {
1156 _loop++;
1157 RecCgemCluster *cluster = (*iter);
1158 int flag=cluster->getflag();
1159 if(flag!=2)continue;
1160
1161 int layerid=cluster->getlayerid();
1162 int sheetid=cluster->getsheetid();
1163 double _phi=cluster->getrecphi();
1164 double _v=cluster->getrecv();
1165 double _Z=cluster->getRecZ();
1166 double _Q=cluster->getenergydeposit();
1167 double t_phi=cluster->getrecphi_mTPC();
1168 double t_v=cluster->getrecv_mTPC();
1169 double t_Z=cluster->getRecZ_mTPC();
1170 int vir_layerid = GetVirLay(layerid, _phi);
1171
1172 if(!(_loop==L_11[i_11]||_loop==L_10[i_10]||_loop==L_00[i_00]||_loop==L_01[i_01]))continue;
1173
1174 Vec_flag.push_back(flag);
1175 Vec_layerid.push_back(layerid);
1176 Vec_Virlayerid.push_back(vir_layerid);
1177 Vec_sheetid.push_back(sheetid);
1178 Vec_layer.push_back(R_layer[layerid]);
1179
1180 if(Switch_CCmTPC==1)
1181 {
1182 Vec_phi.push_back(t_phi);
1183 Vec_v.push_back(t_v);
1184 Vec_Z.push_back(t_Z);
1185
1186 Vec_m_phi.push_back(_phi);
1187 Vec_m_v.push_back(_v);
1188 Vec_m_Z.push_back(_Z);
1189 }
1190 else if(Switch_CCmTPC==0)
1191 {
1192 Vec_phi.push_back(_phi);
1193 Vec_v.push_back(_v );
1194 Vec_Z.push_back(_Z);
1195
1196 Vec_m_phi.push_back(t_phi);
1197 Vec_m_v.push_back(t_v);
1198 Vec_m_Z.push_back(t_Z);
1199
1200 }
1201
1202 }// one combo
1203
1204 if(Vec_layerid.size()!=4&&TEST_N==0)continue;
1205 else if(Vec_layerid.size()!=3&&TEST_N>0)continue;
1206
1207 OrderClusters();
1208
1209 if(Vec_v[0]==Vec_v[1]||Vec_v[0]==Vec_v[2]||Vec_v[0]==Vec_v[3]||Vec_v[1]==Vec_v[2]||Vec_v[1]==Vec_v[3]||Vec_v[2]==Vec_v[3])continue;
1210 StraightLine* l_outer_tmp;
1211
1212 if(TEST_N==1||TEST_N==4) l_outer_tmp=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1213 else if(TEST_N==2||TEST_N==3) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1214 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1215
1216 TMinuit* _fit=Fit(l_outer_tmp,fa,2);
1217 double _chi=_fit->fAmin;
1218 if(_chi<Min_chi)
1219 {
1220 Min_chi=_chi;
1221 C_00=L_00[i_00];
1222 C_01=L_01[i_01];
1223 C_10=L_10[i_10];
1224 C_11=L_11[i_11];
1225 }
1226 delete _fit;
1227 delete l_outer_tmp;
1228 }
1229 }
1230 }
1231 }
1232
1233
1234 Vec_layer.clear();
1235 Vec_phi.clear();
1236 Vec_Z.clear();
1237 Vec_v.clear();
1238 Vec_layerid.clear();
1239 Vec_Virlayerid.clear();
1240 Vec_flag.clear();
1241 Vec_Q.clear();
1242 Vec_sheetid.clear();
1243 Vec_m_layer.clear();
1244 Vec_m_phi.clear();
1245 Vec_m_Z.clear();
1246 Vec_m_v.clear();
1247 Vec_m_layerid.clear();
1248 Vec_m_flag.clear();
1249 Vec_m_Q.clear();
1250 Vec_m_sheetid.clear();
1251 _loop=0;
1252 ic=0;
1253 iter =recCgemClusterCol->begin();
1254 for(;iter!=recCgemClusterCol->end();iter++)
1255 {
1256 _loop++;
1257 RecCgemCluster *cluster = (*iter);
1258 int flag=cluster->getflag();
1259 int layerid=cluster->getlayerid();
1260 int sheetid=cluster->getsheetid();
1261 double _phi=cluster->getrecphi();
1262 double _v=cluster->getrecv();
1263 double _Z=cluster->getRecZ();
1264 double _Q=cluster->getenergydeposit();
1265 double t_phi=cluster->getrecphi_mTPC();
1266 double t_v=cluster->getrecv_mTPC();
1267 double t_Z=cluster->getRecZ_mTPC();
1268 int clusterid=cluster->getclusterid();
1269 int clusterflagb=cluster->getclusterflagb();
1270 int clusterflage=cluster->getclusterflage();
1271 double x=R_layer[layerid]*cos(_phi);
1272 double y=R_layer[layerid]*sin(_phi);
1273 double z=cluster->getRecZ();
1274 int vir_layerid = GetVirLay(layerid, _phi);
1275
1276 if(_loop>MAX_Clusters) break;
1277 if(flag!=2)continue;
1278
1279 tr_id_cluster[ic]=clusterid;
1280 tr_x_cluster[ic]=x;
1281 tr_y_cluster[ic]=y;
1282 tr_z_cluster[ic]=z;
1283 tr_Q_cluster[ic]=_Q;
1284 tr_phi_cluster[ic]=_phi;
1285 tr_layer_cluster[ic]=layerid;
1286 tr_sheet_cluster[ic]=sheetid;
1287 tr_Is_selected_cluster[ic] =0;
1288 ic++;
1289 if(!(_loop==C_00||_loop==C_10||_loop==C_01||_loop==C_11))continue;
1290
1291 tr_Is_selected_cluster[ic-1] =1;
1292
1293 Vec_flag.push_back(flag);
1294 Vec_layerid.push_back(layerid);
1295 Vec_Virlayerid.push_back(vir_layerid);
1296 Vec_sheetid.push_back(sheetid);
1297 Vec_layer.push_back(R_layer[layerid]);
1298 if(Switch_CCmTPC==1)
1299 {
1300 Vec_phi.push_back(t_phi);
1301 Vec_v.push_back(t_v);
1302 Vec_Z.push_back(t_Z);
1303 Vec_m_phi.push_back(_phi);
1304 Vec_m_v.push_back(_v);
1305 Vec_m_Z.push_back(_Z);
1306 }
1307 else if(Switch_CCmTPC==0)
1308 {
1309 Vec_phi.push_back(_phi);
1310 Vec_v.push_back(_v );
1311 Vec_Z.push_back(_Z);
1312
1313 Vec_m_phi.push_back(t_phi);
1314 Vec_m_v.push_back(t_v);
1315 Vec_m_Z.push_back(t_Z);
1316
1317 }
1318
1319 cluster->setTrkId(0);
1320 vecclusters.push_back(cluster);
1321
1322 }// one combo
1323 NC = vecclusters.size();
1324 tr_ncluster = ic;
1325
1326 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
1327 {
1328 if(Vec_layerid[i_cl]==2) _clst_2++;
1329 if(Vec_layerid[i_cl]==1) _clst_1++;
1330 if(Vec_layerid[i_cl]==0) _clst_0++;
1331 }
1332
1333 if(Vec_layerid.size()!=4&&TEST_N==0) return false;
1334 else if(Vec_layerid.size()!=3&&TEST_N>0) return false;
1335
1336 OrderClusters();
1337
1338 if(TEST_N==1||TEST_N==4) l_outer=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1339 else if(TEST_N==2||TEST_N==3) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1340 else if(TEST_N==0) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1341
1342 if(Min_chi<Chisq_cut)
1343 return true;
1344 else return false;
1345}
1346
1347//==============================================================================================================================
1348// Method4: loop the largest N energetic clusters on each hemisphere, then select the combination giving the smallest chisq
1349//==============================================================================================================================
1350
1351bool CgemLineFit::Loop_MaxQ() {
1352 double cl_00(0),cl_01(0),cl_11(0),cl_10(0);
1353 int _loop(0);
1354 double maxQ_11(0),maxQ_10(0),maxQ_00(0);
1355 int max_11(-1),max_10(-1),max_00(-1);
1356 double C_00(0),C_10(0),C_11(0),C_01(0);
1357 vector<double>Vec_00_layer,Vec_00_phi,Vec_00_Z,Vec_00_layerid,Vec_00_v,Vec_00_flag,Vec_00_Q,Vec_00_sheetid;//
1358 vector<double>Vec_01_layer,Vec_01_phi,Vec_01_Z,Vec_01_layerid,Vec_01_v,Vec_01_flag,Vec_01_Q,Vec_01_sheetid;//
1359 vector<double>Vec_10_layer,Vec_10_phi,Vec_10_Z,Vec_10_layerid,Vec_10_v,Vec_10_flag,Vec_10_Q,Vec_10_sheetid;//
1360 vector<double>Vec_11_layer,Vec_11_phi,Vec_11_Z,Vec_11_layerid,Vec_11_v,Vec_11_flag,Vec_11_Q,Vec_11_sheetid;//
1361
1362 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
1363 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
1364 vector<int>L_00,L_01,L_10,L_11;
1365 vector<double>Q_00,Q_01,Q_10,Q_11;
1366 int ic(0);
1367 if(debug) cout<<"start Loop_MaxQ method"<<endl;
1368
1369 for(;iter!=recCgemClusterCol->end();iter++)
1370 {
1371 _loop++;
1372 RecCgemCluster *cluster = (*iter);
1373
1374 int flag=cluster->getflag();
1375 int layerid=cluster->getlayerid();
1376 double _Q=cluster->getenergydeposit();
1377 double _phi=cluster->getrecphi();
1378 double _v=cluster->getrecv();
1379 int sheetid=cluster->getsheetid();
1380 int clusterid=cluster->getclusterid();
1381 int clusterflagb=cluster->getclusterflagb();
1382 int clusterflage=cluster->getclusterflage();
1383 double x=R_layer[layerid]*cos(_phi);
1384 double y=R_layer[layerid]*sin(_phi);
1385 double z=cluster->getRecZ();
1386 if(debug) cout<<"layerid is "<<layerid<<endl;
1387 if(_loop>MAX_Clusters) break;
1388 ic++;
1389 if(flag!=0)continue;
1390
1391 if(layerid==1&&sheetid==1){
1392 L_11.push_back(_loop); Q_11.push_back(_Q); }
1393 if(layerid==1&&sheetid==0){
1394 L_10.push_back(_loop); Q_10.push_back(_Q); }
1395 if(layerid==0&&sheetid==0&&_phi>0){
1396 L_01.push_back(_loop); Q_01.push_back(_Q); }
1397 if(layerid==0&&sheetid==0&&_phi<0){
1398 L_00.push_back(_loop); Q_00.push_back(_Q); }
1399
1400 }
1401
1402 // chose the N X clusters with the largest charge and save corresponding index
1403 vector<int>L_00_s,L_01_s,L_10_s,L_11_s;
1404
1405 L_00_s = GetNMaxQ(Q_00,L_00,Nmax);
1406 L_01_s = GetNMaxQ(Q_01,L_01,Nmax);
1407 L_10_s = GetNMaxQ(Q_10,L_10,Nmax);
1408 L_11_s = GetNMaxQ(Q_11,L_11,Nmax);
1409
1410 if(TEST_N==1){L_11_s.clear();L_11_s.push_back(-1);}
1411 else if(TEST_N==2){L_01_s.clear();L_01_s.push_back(-1);}
1412 else if(TEST_N==3){L_00_s.clear();L_00_s.push_back(-1);}
1413 else if(TEST_N==4){L_10_s.clear();L_10_s.push_back(-1);}
1414
1415 if(debug) cout<<"finish the first loop"<<endl;
1416
1417 double Min_chi(1E10);
1418 NXComb = L_00_s.size()*L_01_s.size()*L_10_s.size()*L_11_s.size();
1419
1420 for(int i_11=0;i_11<L_11_s.size();i_11++){
1421 for(int i_10=0;i_10<L_10_s.size();i_10++){
1422 for(int i_00=0;i_00<L_00_s.size();i_00++){
1423 for(int i_01=0;i_01<L_01_s.size();i_01++){
1424
1425 Vec_layer.clear();
1426 Vec_phi.clear();
1427 Vec_Z.clear();
1428 Vec_v.clear();
1429 Vec_layerid.clear();
1430 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
1431 Vec_flag.clear();
1432 Vec_Q.clear();
1433 Vec_sheetid.clear();
1434 Vec_m_layer.clear();
1435 Vec_m_phi.clear();
1436 Vec_m_Z.clear();
1437 Vec_m_v.clear();
1438 Vec_m_layerid.clear();
1439 Vec_m_flag.clear();
1440 Vec_m_Q.clear();
1441 Vec_m_sheetid.clear();
1442 _loop=0;
1443 iter =recCgemClusterCol->begin();
1444 for(;iter!=recCgemClusterCol->end();iter++)
1445 {
1446 _loop++;
1447 RecCgemCluster *cluster = (*iter);
1448 int flag=cluster->getflag();
1449 if(!(_loop==L_11_s[i_11]||_loop==L_10_s[i_10]||_loop==L_00_s[i_00]||_loop==L_01_s[i_01]))continue;
1450 int layerid=cluster->getlayerid();
1451 int sheetid=cluster->getsheetid();
1452 double _phi=cluster->getrecphi();
1453 double _v=cluster->getrecv();
1454 double _Z=cluster->getRecZ();
1455 double _Q=cluster->getenergydeposit();
1456 double t_phi=cluster->getrecphi_mTPC();
1457 double t_v=cluster->getrecv_mTPC();
1458 double t_Z=cluster->getRecZ_mTPC();
1459
1460 int vir_layerid = GetVirLay(layerid, _phi);
1461
1462 Vec_flag.push_back(flag);
1463 Vec_layerid.push_back(layerid);
1464 Vec_Virlayerid.push_back(vir_layerid);
1465 Vec_sheetid.push_back(sheetid);
1466 Vec_layer.push_back(R_layer[layerid]);
1467 if(Switch_CCmTPC==1)
1468 {
1469 Vec_phi.push_back(t_phi);
1470 Vec_v.push_back(t_v);
1471 Vec_Z.push_back(t_Z);
1472
1473 Vec_m_phi.push_back(_phi);
1474 Vec_m_v.push_back(_v);
1475 Vec_m_Z.push_back(_Z);
1476 }
1477 else if(Switch_CCmTPC==0)
1478 {
1479 Vec_phi.push_back(_phi);
1480 Vec_v.push_back(_v );
1481 Vec_Z.push_back(_Z);
1482
1483 Vec_m_phi.push_back(t_phi);
1484 Vec_m_v.push_back(t_v);
1485 Vec_m_Z.push_back(t_Z);
1486
1487 }
1488 }// one combo
1489 if(Vec_layerid.size()!=4&&TEST_N==0)return false;
1490 else if(Vec_layerid.size()!=3&&TEST_N>0)return false;
1491 if(debug) cout<<"before OrderCluster"<<endl;
1492 OrderClusters();
1493 // 0 1 2 3
1494 // order is layer 1 , top, down, layer 0, top down
1495
1496 StraightLine* l_outer_tmp;
1497
1498 if(debug) cout<<"before IniPar"<<endl;
1499 // lay0top virlay lay0down virlay
1500 if(TEST_N==1||TEST_N==4) l_outer_tmp=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1501 else if(TEST_N==2||TEST_N==3) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1502 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1503
1504 if(debug) cout<<"before fit"<<endl;
1505 TMinuit* _fit=Fit(l_outer_tmp,fa,1);
1506 double _chi=_fit->fAmin;
1507 if(debug) cout<<"finish the fit and the chisq is "<<_chi<<endl;
1508 if(debug) cout<<"----------------------------------------"<<endl;
1509 if(_chi<Min_chi)
1510 {
1511 Min_chi=_chi;
1512 if(TEST_N==0) C_00=Vec_phi[3];
1513 else C_00=-99; //Why set C_00 to -99 Guoaq--2020/11/01
1514 C_01=Vec_phi[2]; //It is because there is only 3 clusters selected if TEST_N!=0
1515 C_10=Vec_phi[1];
1516 C_11=Vec_phi[0];
1517
1518 }
1519 delete _fit;
1520 delete l_outer_tmp;
1521 }
1522 }
1523 }
1524 }
1525
1526 if(debug) cout<<"finish the second loop"<<endl;
1527
1528 iter =recCgemClusterCol->begin();
1529 L_00.clear();L_01.clear();L_10.clear();L_11.clear();
1530 L_00_s.clear();L_01_s.clear();L_10_s.clear();L_11_s.clear();
1531 Q_00.clear();Q_01.clear();Q_10.clear();Q_11.clear();
1532 _loop=0;
1533 int ic_tot(0);
1534 for(;iter!=recCgemClusterCol->end();iter++)
1535 {
1536 _loop++;
1537 RecCgemCluster *cluster = (*iter);
1538 int flag=cluster->getflag();
1539 if(flag!=2)continue;
1540 ic_tot++;
1541 int layerid=cluster->getlayerid();
1542 double _Q=cluster->getenergydeposit();
1543 double _phi=cluster->getrecphi();
1544 int sheetid=cluster->getsheetid();
1545 if(_Q<MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
1546 if(C_00==_phi||C_01==_phi||C_11==_phi||C_10==_phi){
1547 if(layerid==1&&sheetid==1){
1548 L_11.push_back(_loop); Q_11.push_back(_Q); }
1549 if(layerid==1&&sheetid==0){
1550 L_10.push_back(_loop); Q_10.push_back(_Q); }
1551 if(layerid==0&&sheetid==0&&_phi>0){
1552 L_01.push_back(_loop); Q_01.push_back(_Q); }
1553 if(layerid==0&&sheetid==0&&_phi<0){
1554 L_00.push_back(_loop); Q_00.push_back(_Q); }
1555 }
1556 }
1557
1558 if(debug) cout<<"finish the third loop"<<endl;
1559
1560 Min_chi=1E10;
1561
1562 clst_11=L_11.size();
1563 clst_01=L_01.size();
1564 clst_10=L_10.size();
1565 clst_00=L_00.size();
1566
1567 L_00_s = GetNMaxQ(Q_00,L_00,Nmax);
1568 L_01_s = GetNMaxQ(Q_01,L_01,Nmax);
1569 L_10_s = GetNMaxQ(Q_10,L_10,Nmax);
1570 L_11_s = GetNMaxQ(Q_11,L_11,Nmax);
1571
1572 if(TEST_N==1) {L_11_s.clear();L_11_s.push_back(-1);}
1573 else if(TEST_N==2){L_01_s.clear();L_01_s.push_back(-1);}
1574 else if(TEST_N==3){L_00_s.clear();L_00_s.push_back(-1);}
1575 else if(TEST_N==4){L_10_s.clear();L_10_s.push_back(-1);}
1576
1577 NVComb = L_00_s.size()*L_01_s.size()*L_10_s.size()*L_11_s.size();
1578
1579 for(int i_11=0;i_11<L_11_s.size();i_11++){
1580 for(int i_10=0;i_10<L_10_s.size();i_10++){
1581 for(int i_00=0;i_00<L_00_s.size();i_00++){
1582 for(int i_01=0;i_01<L_01_s.size();i_01++){
1583
1584 Vec_layer.clear();
1585 Vec_phi.clear();
1586 Vec_Z.clear();
1587 Vec_v.clear();
1588 Vec_layerid.clear();
1589 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
1590 Vec_flag.clear();
1591 Vec_Q.clear();
1592 Vec_sheetid.clear();
1593 Vec_m_layer.clear();
1594 Vec_m_phi.clear();
1595 Vec_m_Z.clear();
1596 Vec_m_v.clear();
1597 Vec_m_layerid.clear();
1598 Vec_m_flag.clear();
1599 Vec_m_Q.clear();
1600 Vec_m_sheetid.clear();
1601 _loop=0;
1602 iter =recCgemClusterCol->begin();
1603 for(;iter!=recCgemClusterCol->end();iter++)
1604 {
1605 _loop++;
1606 RecCgemCluster *cluster = (*iter);
1607 int flag=cluster->getflag();
1608 if(flag!=2)continue;
1609 if(!(_loop==L_11_s[i_11]||_loop==L_10_s[i_10]||_loop==L_00_s[i_00]||_loop==L_01_s[i_01]))continue;
1610
1611 int layerid=cluster->getlayerid();
1612 int sheetid=cluster->getsheetid();
1613 double _phi=cluster->getrecphi();
1614 double _v=cluster->getrecv();
1615 double _Z=cluster->getRecZ();
1616 double _Q=cluster->getenergydeposit();
1617 double t_phi=cluster->getrecphi_mTPC();
1618 double t_v=cluster->getrecv_mTPC();
1619 double t_Z=cluster->getRecZ_mTPC();
1620 int vir_layerid = GetVirLay(layerid, _phi);
1621
1622 Vec_flag.push_back(flag);
1623 Vec_layerid.push_back(layerid);
1624 Vec_sheetid.push_back(sheetid);
1625 Vec_layer.push_back(R_layer[layerid]);
1626 Vec_Virlayerid.push_back(vir_layerid);
1627
1628
1629 if(Switch_CCmTPC==1)
1630 {
1631 Vec_phi.push_back(t_phi);
1632 Vec_v.push_back(t_v);
1633 Vec_Z.push_back(t_Z);
1634
1635 Vec_m_phi.push_back(_phi);
1636 Vec_m_v.push_back(_v);
1637 Vec_m_Z.push_back(_Z);
1638 }
1639 else if(Switch_CCmTPC==0)
1640 {
1641 Vec_phi.push_back(_phi);
1642 Vec_v.push_back(_v );
1643 Vec_Z.push_back(_Z);
1644
1645 Vec_m_phi.push_back(t_phi);
1646 Vec_m_v.push_back(t_v);
1647 Vec_m_Z.push_back(t_Z);
1648
1649 }
1650
1651 }// one combo
1652 if(Vec_layerid.size()!=4&&TEST_N==0)return false;
1653 else if(Vec_layerid.size()!=3&&TEST_N>0)return false;
1654 OrderClusters();
1655
1656 // discard the combination if V strip is re-used !
1657 if(Vec_v[0]==Vec_v[1]||Vec_v[0]==Vec_v[2]||Vec_v[0]==Vec_v[3]||Vec_v[1]==Vec_v[2]||Vec_v[1]==Vec_v[3]||Vec_v[2]==Vec_v[3])continue;
1658 StraightLine* l_outer_tmp;
1659
1660 if(TEST_N==1||TEST_N==4) l_outer_tmp=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1661 else if(TEST_N==2||TEST_N==3) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1662 else if(TEST_N==0) l_outer_tmp=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1663
1664 TMinuit* _fit=Fit(l_outer_tmp,fa,2);
1665 double _chi=_fit->fAmin;
1666 if(_chi<Min_chi)
1667 {
1668 Min_chi=_chi;
1669 C_00=L_00_s[i_00];
1670 C_01=L_01_s[i_01];
1671 C_10=L_10_s[i_10];
1672 C_11=L_11_s[i_11];
1673
1674 }
1675 delete _fit;
1676 delete l_outer_tmp;
1677 }
1678 }
1679 }
1680 }
1681
1682 Vec_layer.clear();
1683 Vec_phi.clear();
1684 Vec_Z.clear();
1685 Vec_v.clear();
1686 Vec_layerid.clear();
1687 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
1688 Vec_flag.clear();
1689 Vec_Q.clear();
1690 Vec_sheetid.clear();
1691 Vec_m_layer.clear();
1692 Vec_m_phi.clear();
1693 Vec_m_Z.clear();
1694 Vec_m_v.clear();
1695 Vec_m_layerid.clear();
1696 Vec_m_flag.clear();
1697 Vec_m_Q.clear();
1698 Vec_m_sheetid.clear();
1699 if(debug) cout<<"finish the forth loop"<<endl;
1700 _loop=0;
1701 ic=0;
1702 iter =recCgemClusterCol->begin();
1703 for(;iter!=recCgemClusterCol->end();iter++)
1704 {
1705 _loop++;
1706 RecCgemCluster *cluster = (*iter);
1707 int flag=cluster->getflag();
1708 int layerid=cluster->getlayerid();
1709 int sheetid=cluster->getsheetid();
1710 double _phi=cluster->getrecphi();
1711 double _v=cluster->getrecv();
1712 double _Z=cluster->getRecZ();
1713 double _Q=cluster->getenergydeposit();
1714 double t_phi=cluster->getrecphi_mTPC();
1715 double t_v=cluster->getrecv_mTPC();
1716 double t_Z=cluster->getRecZ_mTPC();
1717 int clusterid=cluster->getclusterid();
1718 int clusterflagb=cluster->getclusterflagb();
1719 int clusterflage=cluster->getclusterflage();
1720 double x=R_layer[layerid]*cos(_phi);
1721 double y=R_layer[layerid]*sin(_phi);
1722 double z=cluster->getRecZ();
1723 int vir_layerid = GetVirLay(layerid, _phi);
1724 if(_loop>MAX_Clusters) break;
1725 if(flag!=2)continue;
1726
1727 tr_id_cluster[ic]=clusterid;
1728 tr_x_cluster[ic]=x;
1729 tr_y_cluster[ic]=y;
1730 tr_z_cluster[ic]=z;
1731 tr_Q_cluster[ic]=_Q;
1732 tr_phi_cluster[ic]=_phi;
1733 tr_layer_cluster[ic]=layerid;
1734 tr_sheet_cluster[ic]=sheetid;
1735 tr_Is_selected_cluster[ic] =0;
1736
1737 ic++;
1738
1739 if(!(_loop==C_00||_loop==C_10||_loop==C_01||_loop==C_11)) continue;
1740 tr_Is_selected_cluster[ic-1] =1;
1741
1742 Vec_flag.push_back(flag);
1743 Vec_layerid.push_back(layerid);
1744 Vec_Virlayerid.push_back(vir_layerid);
1745 Vec_sheetid.push_back(sheetid);
1746 Vec_layer.push_back(R_layer[layerid]);
1747 if(Switch_CCmTPC==1)
1748 {
1749 Vec_phi.push_back(t_phi);
1750 Vec_v.push_back(t_v);
1751 Vec_Z.push_back(t_Z);
1752
1753 Vec_m_phi.push_back(_phi);
1754 Vec_m_v.push_back(_v);
1755 Vec_m_Z.push_back(_Z);
1756 }
1757 else if(Switch_CCmTPC==0)
1758 {
1759 Vec_phi.push_back(_phi);
1760 Vec_v.push_back(_v );
1761 Vec_Z.push_back(_Z);
1762
1763 Vec_m_phi.push_back(t_phi);
1764 Vec_m_v.push_back(t_v);
1765 Vec_m_Z.push_back(t_Z);
1766
1767 }
1768
1769 cluster->setTrkId(0);
1770 vecclusters.push_back(cluster);
1771
1772 }// one combo
1773 NC = vecclusters.size();
1774 tr_ncluster = ic;
1775 if(debug) cout<<"finish the fifth loop"<<endl;
1776
1777 for (int i_cl=0;i_cl<Vec_layerid.size();i_cl++)
1778 {
1779 if(Vec_layerid[i_cl]==2) _clst_2++;
1780 if(Vec_layerid[i_cl]==1) _clst_1++;
1781 if(Vec_layerid[i_cl]==0) _clst_0++;
1782 }
1783 if(debug) cout<<"Vec_layerid.size is "<<Vec_layerid.size()<<endl;
1784
1785 if(Vec_layerid.size()!=4&&TEST_N==0)return false;
1786 else if(Vec_layerid.size()!=3&&TEST_N>0)return false;
1787 OrderClusters();
1788
1789 if(TEST_N==1||TEST_N==4) l_outer=IniPar(Vec_phi[1],Vec_Z[1],1,Vec_phi[2],Vec_Z[2],0);
1790 else if(TEST_N==2||TEST_N==3) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1791 else if(TEST_N==0) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1792
1793 if(debug) cout<<"Event "<<event<<", Clu. ID are "<<C_00<<", "<<C_01<<", "<<C_10<<", "<<C_11<<", Loop_MaxQ Min_chi = "<<Min_chi<<endl;
1794
1795 if(Min_chi<Chisq_cut)
1796 return true;
1797 else return false;
1798}
1799
1800
1801
1802//========================================================================================================
1803// useful functions
1804//========================================================================================================
1805
1806StatusCode CgemLineFit ::registerTrack(RecMdcTrackCol*& recMdcTrackCol ,RecMdcHitCol*& recMdcHitCol)
1807{
1808 MsgStream log(msgSvc(), name());
1809 StatusCode sc;
1810 IDataProviderSvc* eventSvc = NULL;
1811 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1812 if (eventSvc) {
1813 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1814 } else {
1815 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1816 return StatusCode::FAILURE;
1817 //return StatusCode::SUCCESS;
1818 }
1819 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
1820
1821
1822 // --- register RecMdcTrack
1823 recMdcTrackCol = new RecMdcTrackCol;
1824 DataObject *aRecMdcTrackCol;
1825 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1826 if(aRecMdcTrackCol != NULL) {
1827 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1828 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1829 }
1830
1831 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
1832 if(sc.isFailure()) {
1833 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
1834 return StatusCode::FAILURE;
1835 }
1836 log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
1837
1838
1839 recMdcHitCol = new RecMdcHitCol;
1840 DataObject *aRecMdcHitCol;
1841 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1842 if(aRecMdcHitCol != NULL) {
1843 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1844 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1845 }
1846
1847 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
1848 if(sc.isFailure()) {
1849 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1850 return StatusCode::FAILURE;
1851 }
1852 log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
1853
1854 return StatusCode::SUCCESS;
1855}//end of stroeTracks
1856
1857
1858
1859void CgemLineFit::fcn( Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1860{
1861 double chisq = 0;
1862 for (int i=0;i<Vec_Virlayerid.size();i++)
1863 {
1864
1865 if(Vec_Virlayerid[i]==5&&bool(sheet_flag&f21))continue;
1866 if(Vec_Virlayerid[i]==4&&bool(sheet_flag&f20))continue;
1867 if(Vec_Virlayerid[i]==3&&bool(sheet_flag&f11))continue;
1868 if(Vec_Virlayerid[i]==2&&bool(sheet_flag&f10))continue;
1869 if(Vec_Virlayerid[i]==1&&bool(sheet_flag&f01))continue;
1870 if(Vec_Virlayerid[i]==0&&bool(sheet_flag&f00))continue;
1871 double dx1=dx(Vec_Virlayerid[i],Vec_layer[i],par[0],par[1],dz_set,tanl_set,Vec_phi[i]);
1872 chisq+=dx1/sigma2;
1873 }
1874 f = chisq;
1875}
1876
1877void CgemLineFit::fcn2( Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1878{
1879 double chisq = 0;
1880
1881 for (int i=0;i<Vec_Virlayerid.size();i++)
1882 {
1883 if(Vec_Virlayerid[i]==5&&bool(sheet_flag&f21))continue;
1884 if(Vec_Virlayerid[i]==4&&bool(sheet_flag&f20))continue;
1885 if(Vec_Virlayerid[i]==3&&bool(sheet_flag&f11))continue;
1886 if(Vec_Virlayerid[i]==2&&bool(sheet_flag&f10))continue;
1887 if(Vec_Virlayerid[i]==1&&bool(sheet_flag&f01))continue;
1888 if(Vec_Virlayerid[i]==0&&bool(sheet_flag&f00))continue;
1889
1890 if(Vec_flag[i]!=2)continue;
1891 double dv1=dV(Vec_Virlayerid[i],Vec_layer[i],dr_set,phi0_set,par[0],par[1],Vec_phi[i],Vec_v[i]);
1892 chisq+=dv1/sigma2;
1893 }
1894 f=chisq;
1895}
1896
1897
1898void CgemLineFit::fcn3( Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1899{
1900
1901 double chisq = 0;
1902 double dv1(0),dx1(0) ;
1903 for (int i=0;i<Vec_Virlayerid.size();i++)
1904 {
1905 if(Vec_flag[i]==2){
1906
1907 if(Vec_Virlayerid[i]==5&&bool(sheet_flag&f21))continue;
1908 if(Vec_Virlayerid[i]==4&&bool(sheet_flag&f20))continue;
1909 if(Vec_Virlayerid[i]==3&&bool(sheet_flag&f11))continue;
1910 if(Vec_Virlayerid[i]==2&&bool(sheet_flag&f10))continue;
1911 if(Vec_Virlayerid[i]==1&&bool(sheet_flag&f01))continue;
1912 if(Vec_Virlayerid[i]==0&&bool(sheet_flag&f00))continue;
1913
1914 dv1=dV(Vec_Virlayerid[i],Vec_layer[i],par[0],par[1],par[2],par[3],Vec_phi[i],Vec_v[i]);
1915 chisq+=(dv1/sigma2);
1916 dx1=dx(Vec_Virlayerid[i],Vec_layer[i],par[0],par[1],par[2],par[3],Vec_phi[i]);
1917 chisq+=dx1/sigma2;
1918 }
1919 }
1920 f=chisq;
1921}
1922
1923
1924double CgemLineFit::dx(double vir_layerid ,double R,double dr,double phi0,double dz,double tanl,double x )
1925{
1926 if(dr>R) return 9999999999999;
1927 StraightLine l1(dr,phi0,dz,tanl);
1928 double phiV_up[2],phiV_down[2];
1929 HepPoint3D Up,Down;
1930 double phiV_up1_old[2],phiV_down1_old[2];
1931 HepPoint3D Up1_old,Down1_old;
1932 bool findinter = false;
1933 if(!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1934 else findinter = Mp->getPointAligned(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1935
1936 if(!findinter) return 9999999999999;
1937
1938 phiV_up[0] = fmod((phiV_up[0]),(2.0*TMath::Pi()));
1939 phiV_down[0] = fmod((phiV_down[0]),(2.0*TMath::Pi()));
1940
1941 if(phiV_up[0]<0) phiV_up[0]+=2.0*TMath::Pi();
1942 if(phiV_down[0]<0) phiV_down[0]+=2.0*TMath::Pi();
1943
1944 double dx1=Min_diff(x,phiV_down[0]);
1945 double dx2=Min_diff(x,phiV_up[0]);
1946 double mdv=min(pow(dx1,2.0),pow(dx2,2.0));
1947 return mdv*R*R;
1948}
1949
1950
1951double CgemLineFit::dV(int vir_layerid,double R ,double dr,double phi0,double z0,double tglam,double x,double V)
1952{
1953
1954 HepPoint3D Up,Down;
1955 StraightLine l1(dr,phi0,z0,tglam);
1956 double phiV_up[2],phiV_down[2];
1957 bool findinter = false;
1958 if(!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1959 else findinter = Mp->getPointAligned(vir_layerid,l1,Up,Down,phiV_up,phiV_down);
1960
1961 if(!findinter) return 9999999999999;
1962
1963 double dx1=Min_diff(x,phiV_down[0]);
1964 double dx2=Min_diff(x,phiV_up[0]);
1965
1966 if(dx1<dx2)return pow(fabs(phiV_down[1]-V),2.0);
1967 else return pow(fabs(phiV_up[1]-V),2.0);
1968
1969}
1970
1971TMinuit* CgemLineFit::Fit(StraightLine* l_outer,int _sheet_flag,int typ )
1972{
1973 // typ == 0 : 2D fit on XY plane and SZ plane indipendently, then use the fit results as initial values to perform 3D fit!
1974 // typ == 1 : 2D fit on XY plane
1975 // typ == 2 : 3D fit on XY and SZ plane
1976
1977 double trkpar[4]={-999,-999,-999,-999};
1978 double trkparErr[4]={-999,-999,-999,-999};
1979
1980 sheet_flag=~_sheet_flag;
1981 dz_set=l_outer->dz();
1982 tanl_set=l_outer->tanl();
1983 Int_t ierflg ,npari,nparx,istat1,istat2,istat3;
1984 Double_t fmin,edm,errdef;
1985 Double_t arglist[10];
1986
1987 arglist[0] = 50;
1988 arglist[1] = 0.01;
1989
1990 if(typ==2){
1991 TMinuit *Min3 = new TMinuit(4);
1992 Min3 -> SetFCN(fcn3);
1993 Min3 -> SetPrintLevel(-1);
1994 Min3 -> mnparm(0, "dr" , l_outer->dr(), 0.1 , -1* R_layer[2], R_layer[2], ierflg);
1995 Min3 -> mnparm(1, "phi0" , l_outer->phi0(), 0.01, -100*acos(-1), 100*acos(-1), ierflg); // Why the range of phi0 is so big -- Guoaq--2020/11/01
1996 Min3 -> mnparm(2, "dz" , l_outer->dz(), 0.1 , -0.5*length, 0.5*length, ierflg);
1997 Min3 -> mnparm(3, "tanL" , l_outer->tanl(), 0.01, -9999, 9999, ierflg);
1998 Min3 -> mnexcm("MIGRAD", arglist, 2, ierflg);
1999
2000 return Min3;
2001 }
2002
2003 if(typ==0){
2004
2005 arglist[0] = 2000000;
2006 arglist[1] = 0.0001;
2007 }
2008 TMinuit *Min = new TMinuit(2);
2009 Min -> SetPrintLevel(-1);
2010 Min -> SetFCN(fcn);
2011 Min -> SetErrorDef(1.0); // For Chisq
2012 Min -> mnparm(0, "dr" ,l_outer->dr(), 0.001, -1*R_layer[2], R_layer[2], ierflg);
2013 Min -> mnparm(1, "phi0" ,l_outer->phi0(), 0.001, -100*acos(-1), 100*acos(-1), ierflg);
2014 arglist[0] = 2000000;
2015 arglist[1] = 0.001;
2016 arglist[0] = 20;
2017 arglist[1] = 0.1;
2018 Min -> mnexcm("MIGRAD", arglist, 2, ierflg);
2019 Min -> mnstat(fmin, edm, errdef, npari, nparx, istat1);
2020 Min -> GetParameter(0, trkpar[0], trkparErr[0]);
2021 Min -> GetParameter(1, trkpar[1], trkparErr[1]);
2022
2023 to_0_2pi(trkpar[1]);
2024 dr_set=trkpar[0];
2025 phi0_set=trkpar[1];
2026 chisq_1=Min->fAmin;
2027 if(typ==1)return Min;
2028 else delete Min;
2029
2030 /////////////////////////Fit to z0 & tanlam////////////
2031 TMinuit *Min2 = new TMinuit(2);
2032 Min2 -> SetFCN(fcn2);
2033 Min2 -> SetPrintLevel(-1);
2034 Min2 -> mnparm(0, "dz" , l_outer->dz(), 0.001 , -0.5*length, 0.5*length, ierflg);
2035 Min2 -> mnparm(1, "tanL" , l_outer->tanl(), 0.001 , -9999 , 9999, ierflg);
2036 Min2 -> mnexcm("MIGRAD", arglist, 2, ierflg);
2037 Min2 -> mnstat(fmin, edm, errdef, npari, nparx, istat2);
2038 for(int i=0; i<2; i++){
2039 Min2 -> GetParameter(i, trkpar[i+2], trkparErr[i+2]);
2040 }
2041 chisq_2=Min2->fAmin;
2042
2043 delete Min2;
2044 //////////////////re-Fit 4 parameters to get correct EDM /////////////
2045
2046
2047 TMinuit *Min3 = new TMinuit(4);
2048 Min3 -> SetFCN(fcn3);
2049 Min3 -> SetPrintLevel(-1);
2050 Min3 -> mnparm(0, "dr" , trkpar[0], 0.001 , -1* R_layer[2], R_layer[2], ierflg);
2051 Min3 -> mnparm(1, "phi0" , trkpar[1], 0.001 , -100*acos(-1), 100*acos(-1), ierflg);
2052 Min3 -> mnparm(2, "dz" , trkpar[2], 0.001 , -0.5*length, 0.5*length, ierflg);
2053 Min3 -> mnparm(3, "tanL" , trkpar[3], 0.001 , -9999, 9999, ierflg);
2054 Min3 -> mnexcm("MIGRAD", arglist, 2, ierflg);
2055
2056 if(typ==0) return Min3;
2057 else delete Min3;
2058}
2059
2060
2061void CgemLineFit::OrderClusters ()
2062{
2063 int noswap(0);
2064 while(noswap!=Vec_flag.size()-1)
2065 {
2066 noswap=0;
2067 for(int i=0;i<Vec_flag.size()-1;i++)
2068 {
2069 if(Vec_flag[i]<Vec_flag[i+1])
2070 { Swap(i,i+1);
2071 }
2072 else if(Vec_flag[i]==Vec_flag[i+1]&&Vec_layerid[i]<Vec_layerid[i+1])
2073 { Swap(i,i+1);
2074 }
2075 else if(Vec_flag[i]==Vec_flag[i+1]&&Vec_layerid[i]==Vec_layerid[i+1]&&(Vec_phi[i]<Vec_phi[i+1]))
2076 {
2077 Swap(i,i+1);}
2078 else {noswap++;
2079 }
2080 }
2081 }
2082}
2083
2084void CgemLineFit::Swap (int i, int j)
2085{
2086
2087 swap(Vec_layer[i],Vec_layer[j]);
2088 swap(Vec_sheetid[i],Vec_sheetid[j]);
2089 swap(Vec_phi[i],Vec_phi[j]);
2090 swap(Vec_Z[i],Vec_Z[j]);
2091 swap(Vec_v[i],Vec_v[j]);
2092 swap(Vec_layerid[i],Vec_layerid[j]);
2093 swap(Vec_Virlayerid[i],Vec_Virlayerid[j]);
2094 swap(Vec_flag[i],Vec_flag[j]);
2095 if(Vec_m_phi.size()==Vec_flag.size()&&(Method==1||Method==2)){
2096 swap(Vec_m_phi[i],Vec_m_phi[j]);
2097 swap(Vec_m_Z[i],Vec_m_Z[j]);
2098 swap(Vec_m_v[i],Vec_m_v[j]);
2099 }
2100}
2101
2102
2103
2104
2105bool CgemLineFit::Erase_outer()
2106{
2107
2108 vector<double>_x,_v;
2109 vector<int>_iter;
2110 int wrong(-1);
2111 for(int i=0;i<3;i++)
2112 {
2113 _x.push_back(Vec_phi[i]);
2114 _v.push_back(Vec_v[i]);
2115 _iter.push_back(i);
2116 }
2117 if(_x[0]==_x[1]&&_v[0]==_v[2])wrong=0;
2118 else if(_x[0]==_x[1]&&_v[1]==_v[2])wrong=1;
2119 else if(_x[0]==_x[2]&&_v[0]==_v[1])wrong=0;
2120 else if(_x[0]==_x[2]&&_v[2]==_v[1])wrong=2;
2121 else if(_x[1]==_x[2]&&_v[1]==_v[0])wrong=1;
2122 else if(_x[1]==_x[2]&&_v[2]==_v[0])wrong=2;
2123 else {
2124 return false;
2125 }
2126 erase(_iter[wrong]);
2127 return true;
2128}
2129
2130void CgemLineFit::erase (int i)
2131{
2132 vector<double>::iterator _iter_v;
2133 _iter_v=Vec_v.begin();
2134 Vec_v.erase( _iter_v+i );
2135
2136 vector<double>::iterator _iter_phi;
2137 _iter_phi=Vec_phi.begin();
2138 Vec_phi.erase( _iter_phi+i );
2139
2140 vector<double>::iterator _iter_Z;
2141 _iter_Z=Vec_Z.begin();
2142 Vec_Z.erase( _iter_Z+i );
2143
2144 vector<double>::iterator _iter_layer;
2145 _iter_layer=Vec_layer.begin();
2146 Vec_layer.erase( _iter_layer+i );
2147
2148 vector<double >::iterator _iter_layerid;
2149 _iter_layerid=Vec_layerid.begin();
2150 Vec_layerid.erase( _iter_layerid+i );
2151
2152 vector<double >::iterator _iter_Virlayerid;
2153 _iter_Virlayerid=Vec_Virlayerid.begin();
2154 Vec_Virlayerid.erase( _iter_Virlayerid+i );
2155
2156 vector<double>::iterator _iter_flag;
2157 _iter_flag=Vec_flag.begin();
2158 Vec_flag.erase( _iter_flag+i );
2159
2160}
2161
2162
2163void CgemLineFit::Discard (int layer)
2164{
2165 for(int i=0;i<Vec_layer.size();i++)
2166 {
2167 if(Vec_layerid[i]==layer) erase(i);
2168 }
2169}
2170
2171StraightLine* CgemLineFit::IniPar(double phi1,double z1,int i ,double phi2,double z2,int j)
2172{
2173 int Gi = int(i/2);
2174 int Gj = int(j/2);
2175
2176 HepPoint3D up(R_layer[Gi]*cos(phi1),R_layer[Gi]*sin(phi1),z1);
2177 HepPoint3D down(R_layer[Gj]*cos(phi2),R_layer[Gj]*sin(phi2),z2);
2178
2179 if(Align_Flag)up=Al->point_invTransform(i ,up);// added by Guoaq-2020/03/13
2180 if(Align_Flag)down=Al->point_invTransform(j ,down);
2181
2182 StraightLine* l=new StraightLine(up,down);
2183 return l;
2184
2185}
2186
2187// I guess this function is designed to remove the mis-combination of the 2D clusters on the layer0
2188// I think the idea is to compare the distance of clusters to the intersections of the given line, and chose the pair of closest one
2189void CgemLineFit::Filter(int layerid, StraightLine* l1)
2190{
2191 if(!Layer_cross(layerid,l1))
2192 {
2193 Discard(layerid);
2194 }
2195
2196 vector <int>number;
2197 vector <double>dst;
2198 vector <double>dst_up;
2199 vector <double>dst_down;
2200 vector <double>v_x;
2201 vector <double>v_v;
2202
2203 if(Run10_Flag)
2204 {
2205 for (int i =0; i<Vec_layer.size();i++)
2206 {
2207 if(Vec_layerid[i]!=layerid)continue;
2208 number.push_back(i);
2209 // obtain the distance between the up/down side of crosssections to the given position (phi, V)
2210 dst_up.push_back(Min_dis_up(layerid,Vec_phi[i],Vec_v[i],l1));
2211 dst_down.push_back(Min_dis_down(layerid,Vec_phi[i],Vec_v[i],l1));
2212
2213 }
2214
2215 for (int j =0; j<number.size();j++)
2216 {
2217 int i=number[j];
2218 }
2219
2220 // sort the clusters according to the up side distance
2221 int Loop;
2222 Loop=0;
2223 while(Loop!=number.size()-1){
2224 Loop=0;
2225 for (int j=0 ;j<number.size()-1;j++)
2226 {
2227 if(dst_up[j]>dst_up[j+1])
2228 {
2229 swap(dst_up[j],dst_up[j+1]);
2230 swap(dst_down[j],dst_down[j+1]);
2231 swap(number[j],number[j+1]);
2232 }
2233 else Loop++;
2234 }
2235 }
2236
2237 // set the index of the nearest cluster
2238 for (int j =0; j<number.size();j++)
2239 {
2240 int i=number[j];
2241 }
2242
2243 // remove the first member
2244 number.erase(number.begin());
2245 dst_up.erase(dst_up.begin());
2246 dst_down.erase(dst_down.begin());
2247
2248 // what's the difference to previous operation?
2249 for (int j =0; j<number.size();j++)
2250 {
2251 int i=number[j];
2252 }
2253
2254 // sort the clusters according to the down side distance
2255 Loop=0;
2256 while(Loop!=number.size()-1){
2257 Loop=0;
2258 for (int j=0 ;j<number.size()-1;j++)
2259 {
2260 if(dst_down[j]>dst_down[j+1])
2261 {
2262 swap(dst_down[j],dst_down[j+1]);
2263 swap(number[j],number[j+1]);
2264 }
2265 else Loop++;
2266 }
2267 }
2268
2269 for (int j =0; j<number.size();j++)
2270 {
2271 int i=number[j];
2272 }
2273
2274 number.erase(number.begin());
2275
2276 }
2277 else // if run10_flag is false
2278 {
2279 for (int i =0; i<Vec_layer.size();i++)
2280 {
2281 if(Vec_layerid[i]!=layerid)continue;
2282 number.push_back(i);
2283
2284 dst.push_back(Min_dis(layerid,Vec_phi[i],Vec_v[i],l1));
2285
2286 }
2287
2288 int Loop(0);
2289
2290 while(Loop!=number.size()-1){
2291 Loop=0;
2292 for (int j=0 ;j<number.size()-1;j++)
2293 {
2294 if(dst[j]>dst[j+1])
2295 {
2296 swap(dst[j],dst[j+1]);
2297 swap(number[j],number[j+1]);
2298 }
2299 else Loop++;
2300 }
2301 }
2302
2303
2304 number.erase(number.begin());
2305
2306 } // end of if run10_flag
2307
2308 for (int j =0; j<number.size();j++)
2309 {
2310 int i=number[j];
2311 }
2312
2313
2314 sort(number.begin(),number.end(),less<int>());
2315 for(int i=number.size()-1;i!=-1;i--)
2316 {
2317 erase(number[i]);
2318 }
2319
2320 return;
2321}
2322
2323
2324
2325
2326bool CgemLineFit::Get_MCInf(){
2327 vector<double>Vec_truth;
2328 Vec_truth.clear();
2329 Vec_truth=MC_truth();
2330
2331 vector<double> _Vec_truth=Trans(Vec_truth[0],Vec_truth[1],Vec_truth[2],Vec_truth[3]);
2332
2333 MC_par0=_Vec_truth[0];
2334 MC_par1=_Vec_truth[1];
2335 MC_par2=_Vec_truth[2];
2336 MC_par3=_Vec_truth[3];
2337 MC_px=Vec_truth[4];
2338 MC_py=Vec_truth[5];
2339 MC_pz=Vec_truth[6];
2340
2341 vector<double>V_MC_clusters=Get_Clusters(Vec_truth);
2342 x_0_up_mc=V_MC_clusters[0];
2343 v_0_up_mc=V_MC_clusters[1];
2344 x_0_down_mc=V_MC_clusters[2];
2345 v_0_down_mc=V_MC_clusters[3];
2346
2347 x_1_up_mc=V_MC_clusters[4];
2348 v_1_up_mc=V_MC_clusters[5];
2349 x_1_down_mc=V_MC_clusters[6];
2350 v_1_down_mc=V_MC_clusters[7];
2351
2352 x_2_up_mc=V_MC_clusters[8];
2353 v_2_up_mc=V_MC_clusters[9];
2354 x_2_down_mc=V_MC_clusters[10];
2355 v_2_down_mc=V_MC_clusters[11];
2356 return true;
2357}
2358
2359
2360vector<double>CgemLineFit::MC_truth()
2361{
2362 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
2363
2364 Event::McParticleCol::iterator iter_mc=mcParticleCol->begin();
2365 HepLorentzVector p4_mu(0,0,0,0);
2366 HepLorentzVector p4_pos(0,0,0,0);
2367 for(;iter_mc!=mcParticleCol->end();iter_mc++)
2368 {
2369 if (!(*iter_mc)->decayFromGenerator()) continue;
2370 HepLorentzVector p4_mc(0,0,0,0);
2371
2372 if (13==(*iter_mc)->particleProperty())
2373 {
2374
2375 p4_mu = (*iter_mc)->initialFourMomentum();
2376 p4_pos = (*iter_mc)->initialPosition();
2377 }
2378 }
2379 vector<double>_mc;
2380 Hep3Vector _BP(p4_pos.x()*10,p4_pos.y()*10,p4_pos.z()*10);
2381 pos_x=_BP[0];
2382 pos_y=_BP[1];
2383 pos_z=_BP[2];
2384 _mc=Get4Par(p4_mu,_BP);
2385
2386 _mc.push_back(p4_mu.px());
2387 _mc.push_back(p4_mu.py());
2388 _mc.push_back(p4_mu.pz());
2389 if(HitPos(p4_mu,_BP))
2390 _mc.push_back(1);
2391 else
2392 _mc.push_back(0);
2393
2394 return _mc;
2395}
2396
2397
2398
2399
2400vector<double>CgemLineFit::Get4Par(HepLorentzVector p4,Hep3Vector bp)
2401{
2402 vector<double>_V_par;
2403 _V_par.clear();
2404 double xb=p4.px();
2405 double yb=p4.py();
2406 double zb=p4.pz();
2407 double xa=bp[0];
2408 double ya=bp[1];
2409 double za=bp[2];
2410 double k=-1*(xa*xb+ya*yb)/(xb*xb+yb*yb);
2411 double xc=k*xb+xa;
2412 double yc=k*yb+ya;
2413 double zc=k*zb+za;
2414 HepPoint3D a(xa,ya,za);
2415 HepPoint3D c(xc,yc,zc);
2416 StraightLine l(a,c);
2417 _V_par.push_back(l.dr());
2418 _V_par.push_back(l.phi0());
2419 _V_par.push_back(l.dz());
2420 _V_par.push_back(l.tanl());
2421 return _V_par;
2422}
2423
2424
2425bool CgemLineFit::HitPos(HepLorentzVector p4,Hep3Vector bp)
2426{
2427 double L(500),W(160),H(600);
2428
2429 vector<double>_V_par;
2430 _V_par.clear();
2431 double xb=p4.px();
2432 double yb=p4.py();
2433 double zb=p4.pz();
2434 double xa=bp[0];
2435 double ya=bp[1];
2436 double za=bp[2];
2437 double k=H/yb*-1;
2438 double xc=k*xb+xa;
2439 double yc=k*yb+ya;
2440 double zc=k*zb+za;
2441 hit_x=xc;
2442 hit_y=yc;
2443 hit_z=zc;
2444
2445 if(fabs(xc)>W/2)return false;
2446 if(fabs(zc)>L/2)return false;
2447
2448 else return true;
2449}
2450
2451
2452
2453void CgemLineFit::Get_OtherIniPar(int clst_0, int clst_1, int clst_2)
2454{
2455
2456 if(clst_2==0)
2457 {
2458
2459 if(clst_0>=2&&clst_1>=2){
2460 StraightLine* l_ini_i=IniPar(Vec_phi[2],Vec_Z[2],1,Vec_phi[3],Vec_Z[3],0); // why use index of 4 and 5 because we only have 2 layer. I guess hongpeng want to get the initial parameters by the innermost layer
2461
2462 vector<double> inii_par=Trans(l_ini_i->dr(),l_ini_i->phi0(),l_ini_i->dz(),l_ini_i->tanl());
2463 inii_0=inii_par[0];
2464 inii_1=inii_par[1];
2465 inii_2=inii_par[2];
2466 inii_3=inii_par[3];
2467 delete l_ini_i;
2468 }
2469 if(clst_1>=2){
2470 StraightLine* l_ini_m=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2); // I guess Hongpeng want to get the initial parameters by the middle layer
2471 vector<double> inim_par=Trans(l_ini_m->dr(),l_ini_m->phi0(),l_ini_m->dz(),l_ini_m->tanl());
2472
2473 inim_0=inim_par[0];
2474 inim_1=inim_par[1];
2475 inim_2=inim_par[2];
2476 inim_3=inim_par[3];
2477 delete l_ini_m;
2478 }
2479
2480 }
2481 else
2482 {
2483 if(clst_0>=2&&clst_1>=2){
2484 StraightLine* l_ini_i=IniPar(Vec_phi[4],Vec_Z[4],1,Vec_phi[5],Vec_Z[5],0); // why use index of 4 and 5 because we only have 2 layer. I guess hongpeng want to get the initial parameters by the innermost layer
2485
2486 vector<double> inii_par=Trans(l_ini_i->dr(),l_ini_i->phi0(),l_ini_i->dz(),l_ini_i->tanl());
2487 inii_0=inii_par[0];
2488 inii_1=inii_par[1];
2489 inii_2=inii_par[2];
2490 inii_3=inii_par[3];
2491 delete l_ini_i;
2492 }
2493 if(clst_1>=2){
2494 StraightLine* l_ini_m=IniPar(Vec_phi[2],Vec_Z[2],3,Vec_phi[3],Vec_Z[3],2); // I guess Hongpeng want to get the initial parameters by the middle layer
2495 vector<double> inim_par=Trans(l_ini_m->dr(),l_ini_m->phi0(),l_ini_m->dz(),l_ini_m->tanl());
2496
2497 inim_0=inim_par[0];
2498 inim_1=inim_par[1];
2499 inim_2=inim_par[2];
2500 inim_3=inim_par[3];
2501 delete l_ini_m;
2502 }
2503 }
2504}
2505
2506
2507void CgemLineFit::Rec_Clusters()
2508{
2509
2510 for (int i=0;i<Vec_Virlayerid.size();i++)
2511 {
2512
2513 if(Vec_Virlayerid[i]==5)
2514 {
2515 x_2_up = Vec_phi[i];
2516 v_2_up = Vec_v[i];
2517 z_2_up = Vec_Z[i];
2518 }
2519 if(Vec_Virlayerid[i]==4)
2520 {
2521 x_2_down = Vec_phi[i];
2522 v_2_down = Vec_v[i];
2523 z_2_down = Vec_Z[i];
2524 }
2525 if(Vec_Virlayerid[i]==3)
2526 {
2527 x_1_up = Vec_phi[i];
2528 v_1_up = Vec_v[i];
2529 z_1_up = Vec_Z[i];
2530 }
2531
2532 if(Vec_Virlayerid[i]==2)
2533 {
2534 x_1_down = Vec_phi[i];
2535 v_1_down = Vec_v[i];
2536 z_1_down = Vec_Z[i];
2537 }
2538
2539 if(Vec_Virlayerid[i]==1)
2540 {
2541 x_0_up = Vec_phi[i];
2542 v_0_up = Vec_v[i];
2543 z_0_up = Vec_Z[i];
2544 }
2545
2546 if(Vec_Virlayerid[i]==0)
2547 {
2548 x_0_up = Vec_phi[i];
2549 v_0_up = Vec_v[i];
2550 z_0_up = Vec_Z[i];
2551 }
2552
2553 }
2554
2555}
2556
2557
2558
2559void CgemLineFit::Rec_Clusters_mTPC()
2560{
2561
2562 for (int i=0;i<Vec_Virlayerid.size();i++)
2563 {
2564
2565 if(Vec_Virlayerid[i]==5)
2566 {
2567 x_2_up_m = Vec_m_phi[i];
2568 v_2_up_m = Vec_m_v[i];
2569 z_2_up_m = Vec_m_Z[i];
2570 }
2571 if(Vec_Virlayerid[i]==4)
2572 {
2573 x_2_down_m = Vec_m_phi[i];
2574 v_2_down_m = Vec_m_v[i];
2575 z_2_down_m = Vec_m_Z[i];
2576 }
2577 if(Vec_Virlayerid[i]==3)
2578 {
2579 x_1_up_m = Vec_m_phi[i];
2580 v_1_up_m = Vec_m_v[i];
2581 z_1_up_m = Vec_m_Z[i];
2582 }
2583
2584 if(Vec_Virlayerid[i]==2)
2585 {
2586 x_1_down_m = Vec_m_phi[i];
2587 v_1_down_m = Vec_m_v[i];
2588 z_1_down_m = Vec_m_Z[i];
2589 }
2590
2591 if(Vec_Virlayerid[i]==1)
2592 {
2593 x_0_up_m = Vec_m_phi[i];
2594 v_0_up_m = Vec_m_v[i];
2595 z_0_up_m = Vec_m_Z[i];
2596 }
2597
2598 if(Vec_Virlayerid[i]==0)
2599 {
2600 x_0_up_m = Vec_m_phi[i];
2601 v_0_up_m = Vec_m_v[i];
2602 z_0_up_m = Vec_m_Z[i];
2603 }
2604
2605 }
2606
2607}
2608
2609
2610void CgemLineFit::Fit_Clusters(double par[])
2611{
2612 vector<double>Vec_truth;
2613 Vec_truth.push_back(par[0]);
2614 Vec_truth.push_back(par[1]);
2615 Vec_truth.push_back(par[2]);
2616 Vec_truth.push_back(par[3]);
2617
2618 vector<double>clusters=Get_Clusters(Vec_truth);
2619
2620 x_0_up_f= clusters[0];
2621 v_0_up_f= clusters[1];
2622 x_0_down_f= clusters[2];
2623 v_0_down_f= clusters[3];
2624
2625 x_1_up_f= clusters[4];
2626 v_1_up_f= clusters[5];
2627 x_1_down_f= clusters[6];
2628 v_1_down_f= clusters[7];
2629
2630 x_2_up_f= clusters[8];
2631 v_2_up_f= clusters[9];
2632 x_2_down_f= clusters[10];
2633 v_2_down_f= clusters[11];
2634
2635}
2636
2637void CgemLineFit::GetIntersection()
2638{
2639
2640 if(debug) cout<<"Start Getintersection"<<endl;
2641 Double_t fmin,edm,errdef;
2642 Int_t ierflg ,npari,nparx,istat;
2643 double trkpar[4]={-999,-999,-999,-999};
2644 double trkparErr[4]={-999,-999,-999,-999};
2645
2646 double phiV_up[2],phiV_down[2];
2647 HepPoint3D Up,Down;
2648 double phiV_up_Tmp[2],phiV_down_Tmp[2];
2649 HepPoint3D Up_Tmp,Down_Tmp;
2650
2651 double Ang[3];
2652 StraightLine* l1a=IniPar(Vec_phi[2],Vec_Z[2],1,Vec_phi[3],Vec_Z[3],0);
2653 TMinuit* fit=Fit(l1a,fa-f11,0);
2654 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat);
2655 inter_1_flag=istat;
2656 for(int i=0; i<4; i++){
2657 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2658 StraightLine* l1b=new StraightLine(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
2659 Store_Trk(fit,1000);
2660 delete fit;
2661
2662 if(Align_FLAG)
2663 {
2664 Mp->getPointAligned(2,l1b,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2665 Mp->getPointAligned(3,l1b,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2666 }
2667 else Mp->getPointIdealGeom(1*2,l1b,Up,Down,phiV_up,phiV_down);
2668
2669 to_0_2pi(phiV_up[0]);
2670 to_0_2pi(phiV_down[0]);
2671 inter_1_x=phiV_up[0];
2672 inter_1_V=phiV_up[1];
2673 inter_1_z=Up[2];
2674
2675
2676 inter_4_x=phiV_down[0];
2677 inter_4_V=phiV_down[1];
2678 inter_4_z=Down[2];
2679
2680
2681 HepPoint3D l1_up(cos(Vec_phi[0])*R_layer[1],sin(Vec_phi[0])*R_layer[1],Vec_Z[0]);
2682 HepPoint3D l1_down(cos(Vec_phi[1])*R_layer[1],sin(Vec_phi[1])*R_layer[1],Vec_Z[1]);
2683
2684 IncidentAngle(l1b,l1_up,pl_10->getStereoAngle(),Ang);
2685 ang_1=Ang[0];
2686 ang_1_x=Ang[1];
2687 ang_1_v=Ang[2];
2688
2689 IncidentAngle(l1b,l1_down,pl_11->getStereoAngle(),Ang);
2690 ang_4=Ang[0];
2691 ang_4_x=Ang[1];
2692 ang_4_v=Ang[2];
2693
2694 delete l1a;
2695 delete l1b;
2696
2697 StraightLine * l2a=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
2698 fit=Fit(l2a,fa-f01,0);
2699 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat);
2700 inter_2_flag=istat;
2701 for(int i=0; i<4; i++){
2702 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2703 Store_Trk(fit,2000);
2704 delete fit;
2705 StraightLine* l2b=new StraightLine(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
2706 if(Align_FLAG)
2707 {
2708 Mp->getPointAligned(0,l2b,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2709 Mp->getPointAligned(1,l2b,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2710 }
2711 else Mp->getPointIdealGeom(0*2,l2b,Up,Down,phiV_up,phiV_down);
2712 to_0_2pi(phiV_up[0]);
2713 inter_2_x=phiV_up[0];
2714 inter_2_V=phiV_up[1];
2715 inter_2_z=Up[2];
2716 HepPoint3D l0_up(cos(Vec_phi[2])*R_layer[0],sin(Vec_phi[2])*R_layer[0],Vec_Z[2]);
2717 IncidentAngle(l2b,l0_up,pl_00->getStereoAngle(),Ang);
2718 ang_2=Ang[0];
2719 ang_2_x=Ang[1];
2720 ang_2_v=Ang[2];
2721 delete l2a;
2722 delete l2b;
2723
2724 if(Run10_Flag)
2725 {
2726
2727 StraightLine * l3a=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
2728 fit=Fit(l3a,fa-f00,0);
2729 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat);
2730 inter_3_flag=istat;
2731 for(int i=0; i<4; i++){
2732 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2733 Store_Trk(fit,3000);
2734 delete fit;
2735 StraightLine* l3b=new StraightLine(trkpar[0],trkpar[1],trkpar[2],trkpar[3]);
2736
2737 if(Align_FLAG)
2738 {
2739 Mp->getPointAligned(0,l3b,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2740 Mp->getPointAligned(1,l3b,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2741 }
2742 else Mp->getPointIdealGeom(0*2,l3b,Up,Down,phiV_up,phiV_down);
2743 to_0_2pi(phiV_down[0]);
2744 inter_3_x=phiV_down[0];
2745 inter_3_V=phiV_down[1];
2746 inter_3_z=Down[2];
2747
2748 HepPoint3D l0_down(cos(Vec_phi[1])*R_layer[1],sin(Vec_phi[1])*R_layer[1],Vec_Z[1]);
2749 IncidentAngle(l3b,l0_down,pl_00->getStereoAngle(),Ang);
2750 ang_3=Ang[0];
2751 ang_3_x=Ang[1];
2752 ang_3_v=Ang[2];
2753 delete l3a;
2754 delete l3b;
2755 }
2756
2757}
2758
2759vector<double>CgemLineFit::Get_Clusters(vector<double>Vec_truth)
2760{
2761 vector<double>clusters;
2762 StraightLine l(Vec_truth[0],Vec_truth[1],Vec_truth[2],Vec_truth[3]);
2763 double phiV_up[2],phiV_down[2];
2764 HepPoint3D Up,Down;
2765
2766 double phiV_up_Tmp[2],phiV_down_Tmp[2];
2767 HepPoint3D Up_Tmp,Down_Tmp;
2768
2769 if(!Align_Flag) Mp->getPointIdealGeom(0*2,l,Up,Down,phiV_up,phiV_down);
2770 else
2771 {
2772 Mp->getPointAligned(0,l,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2773 Mp->getPointAligned(1,l,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2774 }
2775 to_0_2pi(phiV_up[0]);
2776 to_0_2pi(phiV_down[0]);
2777
2778 clusters.push_back(phiV_up[0]);
2779 clusters.push_back(phiV_up[1]);
2780 clusters.push_back(phiV_down[0]);
2781 clusters.push_back(phiV_down[1]);
2782
2783 if(!Align_Flag) Mp->getPointIdealGeom(1*2,l,Up,Down,phiV_up,phiV_down);
2784 else
2785 {
2786 Mp->getPointAligned(2,l,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2787 Mp->getPointAligned(3,l,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2788 }
2789 to_0_2pi(phiV_up[0]);
2790 to_0_2pi(phiV_down[0]);
2791
2792 clusters.push_back(phiV_up[0]);
2793 clusters.push_back(phiV_up[1]);
2794 clusters.push_back(phiV_down[0]);
2795 clusters.push_back(phiV_down[1]);
2796
2797 if(!Align_Flag) Mp->getPointIdealGeom(2*2,l,Up,Down,phiV_up,phiV_down);
2798 else
2799 {
2800 Mp->getPointAligned(4,l,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
2801 Mp->getPointAligned(5,l,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
2802 }
2803 to_0_2pi(phiV_up[0]);
2804 to_0_2pi(phiV_down[0]);
2805
2806 clusters.push_back(phiV_up[0]);
2807 clusters.push_back(phiV_up[1]);
2808 clusters.push_back(phiV_down[0]);
2809 clusters.push_back(phiV_down[1]);
2810
2811 return clusters;
2812}
2813
2814void CgemLineFit::to_0_2pi(double &arg)
2815{
2816
2817 arg = fmod((arg),(2.0*TMath::Pi()));
2818 if(arg<0) arg += 2.0*TMath::Pi();
2819}
2820
2821double CgemLineFit::Min_diff(double arg1,double arg2)
2822{
2823 to_0_2pi(arg1);to_0_2pi(arg2);
2824 double diff_raw=fabs(arg1-arg2);
2825
2826 if(diff_raw>acos(-1)) return (2*acos(-1)-diff_raw);
2827 else return diff_raw;
2828}
2829
2830double CgemLineFit::Min_dis(int R_id,double x,double v,StraightLine* l1 )
2831{
2832 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2833 double phiV_up[2],phiV_down[2];
2834 HepPoint3D Up,Down;
2835
2836 int vir_R_id = GetVirLay(R_id, x);
2837
2838 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2839 else Mp->getPointAligned(vir_R_id,l1,Up,Down,phiV_up,phiV_down);
2840
2841 double ds1=
2842 pow((phiV_up[1]-v),2)+
2843 pow(R_layer[R_id]*Min_diff(phiV_up[0],x),2);
2844
2845 double ds2=
2846 pow((phiV_down[1]-v),2)+
2847 pow(R_layer[R_id]*Min_diff(phiV_down[0],x),2);
2848
2849 return min(ds1,ds2);
2850}
2851
2852
2853double CgemLineFit::Min_dis_up(int R_id,double x,double z,StraightLine* l1 )
2854{
2855 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2856 double phiV_up[2],phiV_down[2];
2857 HepPoint3D Up,Down;
2858 int vir_R_id = GetVirLay(R_id, x);
2859 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2860 else Mp->getPointAligned(vir_R_id,l1,Up,Down,phiV_up,phiV_down);
2861 double ds1=
2862 pow((phiV_up[1]-z),2)+
2863 pow(R_layer[R_id]*Min_diff(phiV_up[0],x),2);
2864
2865
2866 return ds1;
2867}
2868
2869
2870double CgemLineFit::Min_dis_down(int R_id,double x,double v,StraightLine* l1 )
2871{
2872
2873 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2874 double phiV_up[2],phiV_down[2];
2875 HepPoint3D Up,Down;
2876 int vir_R_id = GetVirLay(R_id, x);
2877 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2878 else Mp->getPointAligned(vir_R_id,l1,Up,Down,phiV_up,phiV_down);
2879
2880 double ds2=
2881 pow((phiV_down[1]-v),2)+
2882 pow(R_layer[R_id]*Min_diff(phiV_down[0],x),2);
2883
2884 return ds2;
2885}
2886
2887
2888bool CgemLineFit::Layer_cross(int R_id,StraightLine* l1 )
2889{
2890
2891 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
2892 double phiV_up[2],phiV_down[2];
2893
2894 HepPoint3D Up,Down;
2895 if(!Align_Flag) Mp->getPointIdealGeom(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2896 else Mp->getPointAligned(R_id*2,l1,Up,Down,phiV_up,phiV_down);
2897 if(Up[0]>9999||Up[1]>9999||Down[0]>9999||Down[1]>9999)
2898 return false;
2899 else return true;
2900
2901}
2902
2903
2904vector<double> CgemLineFit::Trans(double par0,double par1,double par2,double par3)
2905{
2906 vector<double> rotat;
2907 rotat.clear();
2908
2909 rotat.push_back(par0);
2910 rotat.push_back(par1);
2911 rotat.push_back(par2);
2912 rotat.push_back(par3);
2913 return rotat;
2914
2915}
2916
2917int CgemLineFit::GetVirLay(int geolay, double phi)
2918{
2919 int virsheet = (phi<0)? 0:1;
2920 int virlay = geolay*2 +virsheet;
2921 return virlay;
2922}
2923
2924vector<int> CgemLineFit::GetNMaxQ(vector<double> Q_11,vector<int> L_11,int Nmax=3)
2925{
2926 // chose the N X clusters with the largest charge and save corresponding index
2927 vector<int> L_11_s;
2928 vector<double> Q_11_s;
2929 double t, tt;
2930 if (L_11.size()>Nmax)
2931 {
2932 for(int i=0;i<Nmax;i++)
2933 {
2934 for(int j=i+1;j<L_11.size();j++)
2935 {
2936 if(Q_11[i]<Q_11[j])
2937 {
2938 t=Q_11[i];
2939 Q_11[i]=Q_11[j];
2940 Q_11[j]=t;
2941
2942 tt=L_11[i];
2943 L_11[i]=L_11[j];
2944 L_11[j]=tt;
2945 }
2946 }
2947 Q_11_s.push_back(Q_11[i]);
2948 L_11_s.push_back(L_11[i]);
2949 }
2950
2951 }
2952 else
2953 {
2954 L_11_s.assign(L_11.begin(), L_11.end());
2955 Q_11_s.assign(Q_11.begin(), Q_11.end());
2956 }
2957 return L_11_s;
2958}
2959
2960
2961
2962
2963
2964void CgemLineFit::FilterClusters ()
2965{
2966 int max_1_0(0),max_0_0(0);
2967 int _size=Vec_flag.size();
2968 max_1_0=_size;max_0_0=_size;
2969
2970 for(int i=0;i<Vec_flag.size();i++)
2971 {
2972 if(1==Vec_layerid[i]&&0==Vec_sheetid[i])
2973 {max_1_0==i;break;}
2974 }
2975
2976 for(int j=0;j<Vec_flag.size();j++)
2977 {
2978 if(0==Vec_layerid[j]&&0==Vec_sheetid[j])
2979 {max_0_0==j;break;}
2980 }
2981 for(int k=Vec_flag.size()-1;k>max_0_0;k--)
2982 {erase(k);}
2983 for(int m=Vec_flag.size()-2;m>max_1_0;m--)
2984 {erase(m);}
2985}
2986
2987void CgemLineFit::Store_Trk (TMinuit*fit ,int trkid )
2988{
2989 double trkpar[4]={-999,-999,-999,-999};
2990 double trkparErr[4]={-999,-999,-999,-999};
2991
2992 Int_t ierflg ,npari,nparx,istat3;
2993 Double_t fmin,edm,errdef;
2994
2995 fit -> mnstat(fmin, edm, errdef, npari, nparx, istat3);
2996 for(int i=0; i<4; i++){
2997 fit -> GetParameter(i, trkpar[i], trkparErr[i]);}
2998 double emat[16];
2999 fit->mnemat(emat,4);
3000
3001 RecMdcTrack* recMdcTrack = new RecMdcTrack();
3002
3003 double helixPar[5];
3004 helixPar[0]=trkpar[0];
3005 helixPar[1]=trkpar[1];
3006 helixPar[2]=0.0;
3007 helixPar[3]=trkpar[2];
3008 helixPar[4]=trkpar[3];
3009 double errorMat[15];
3010 int k = 0;
3011 errorMat[0]=emat[0];
3012 errorMat[1]=emat[1];
3013 errorMat[2]=0;
3014 errorMat[3]=emat[2];
3015 errorMat[4]=emat[3];
3016 errorMat[5]=emat[5];
3017 errorMat[6]=0;
3018 errorMat[7]=emat[6];
3019 errorMat[8]=emat[7];
3020 errorMat[9]=0;
3021 errorMat[10]=0;
3022 errorMat[11]=0;
3023 errorMat[12]=emat[10];
3024 errorMat[13]=emat[11];
3025 errorMat[14]=emat[15];
3026
3027 recMdcTrack->setChi2(fit->fAmin);
3028 recMdcTrack->setError(errorMat);
3029 recMdcTrack->setHelix(helixPar);
3030 recMdcTrack->setTrackId(trkid);
3031 if(trkid==0){
3032 recMdcTrack->setNcluster(NC);
3033 recMdcTrack->setVecClusters(vecclusters);
3034 m_trackList_tds->push_back(recMdcTrack);
3035 }
3036
3037}
3038
3039void CgemLineFit::IncidentAngle(StraightLine* cosmic_ray,HepPoint3D cluster ,double theta,double ang[] )
3040{
3041
3042 const Hep3Vector Normals(cluster[0],cluster[1],0);
3043
3044 HepPoint3D x1=cosmic_ray->x(1);
3045 HepPoint3D x2=cosmic_ray->x(-1);
3046 const Hep3Vector CosmicRay(x2[0]-x1[0],x2[1]-x1[1],x2[2]-x1[2]);
3047 const Hep3Vector z(0,0,1);
3048 const Hep3Vector phi=z.cross(Normals);
3049 Hep3Vector _phi=phi;
3050
3051 const Hep3Vector V=_phi.rotate(-1*theta,Normals);
3052
3053 const Hep3Vector N_V=Normals.cross(V);
3054 const Hep3Vector N_phi=Normals.cross(phi);
3055
3056 Hep3Vector C_phi=CosmicRay-CosmicRay.project(N_phi);
3057 ang[0]=Normals.angle(C_phi);
3058
3059 Hep3Vector C_V=CosmicRay-CosmicRay.project(N_V);
3060 ang[1]=Normals.angle(C_V);
3061 ang[2]=Normals.angle(CosmicRay);
3062
3063
3064}
3065
3066
3067
3068
3069
3070void CgemLineFit::GetRes( StraightLine* l,int i_inter)
3071{
3072 int i_layer=0;
3073 int i_sheet=0;
3074 if(i_inter==1||i_inter==4)
3075 {
3076 i_layer=1;
3077 if(i_inter==1) i_sheet=1;
3078 }
3079 double phiV_up[2],phiV_down[2];
3080 HepPoint3D Up,Down;
3081 double phiV_up_Tmp[2],phiV_down_Tmp[2];
3082 HepPoint3D Up_Tmp,Down_Tmp;
3083
3084 if(!Align_Flag) Mp->getPointIdealGeom(i_layer*2,l,Up,Down,phiV_up,phiV_down);
3085 else
3086 {
3087 Mp->getPointAligned(i_layer*2,l,Up_Tmp,Down,phiV_up_Tmp,phiV_down);
3088 Mp->getPointAligned(i_layer*2+1,l,Up,Down_Tmp,phiV_up,phiV_down_Tmp);
3089 }
3090 double V,phi,min_dst(99999),min_V(99999),min_phi(99999);
3091
3092 V=phiV_up[1];
3093 phi=phiV_up[0];
3094
3095 if(i_inter==3||i_inter==4){
3096 V=phiV_down[1];
3097 phi=phiV_down[0];
3098 }
3099
3100 int _loop=0;
3101 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
3102 RecCgemClusterCol::iterator iter =recCgemClusterCol->begin();
3103 int ind=0;
3104 for(;iter!=recCgemClusterCol->end();iter++)
3105 {
3106 _loop++;
3107 RecCgemCluster *cluster = (*iter);
3108 int flag=cluster->getflag();
3109
3110 int layerid=cluster->getlayerid();
3111 if(layerid!=i_layer)continue;
3112 int sheetid=cluster->getsheetid();
3113 if(sheetid!=i_sheet)continue;
3114 double _Q=cluster->getenergydeposit();
3115 double _phi=cluster->getrecphi();
3116 double _V=cluster->getrecv();
3117 double _size=fabs(cluster->getclusterflagb()-cluster->getclusterflage());
3118 if(flag!=2)continue;
3119 double dst=pow(_V-V,2.0)+pow(_phi-phi,2.0);
3120 if(dst<min_dst)
3121 {
3122 min_dst=dst;
3123 min_V=_V;
3124 min_phi=_phi;
3125 }
3126 }
3127
3128 Test_phi=phi;
3129 Test_V=V;
3130 Res_phi=min_phi;
3131 Res_V=min_V;
3132
3133 double Ang[3];
3134
3135}
3136
3137
3138double CgemLineFit::RealPhi(double SimuPhi)
3139{
3140 double RPhi;
3141 if(SimuPhi <= TMath::Pi()) RPhi = SimuPhi;
3142 else RPhi = SimuPhi-2*TMath::Pi();
3143 return RPhi;
3144
3145}
3146
3147
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
Al
CgemGeoAlign * Al
Definition: CgemLineFit.cxx:89
Vec_m_layer
vector< double > Vec_m_layer
Definition: CgemLineFit.cxx:86
Vec_Virlayerid
vector< double > Vec_Virlayerid
Definition: CgemLineFit.cxx:85
f21
int f21(1)
_clst_2
int _clst_2(0)
phi0_set
double phi0_set
Definition: CgemLineFit.cxx:106
GeoLayer2
CgemGeoLayer * GeoLayer2
Definition: CgemLineFit.cxx:92
Vec_Q
vector< double > Vec_Q
Definition: CgemLineFit.cxx:85
f11
int f11(2)
Vec_flag
vector< double > Vec_flag
Definition: CgemLineFit.cxx:85
sigma2
double sigma2(0)
f01
int f01(8)
R_layer
double R_layer[3]
Definition: CgemLineFit.cxx:107
Vec_m_layerid
vector< double > Vec_m_layerid
Definition: CgemLineFit.cxx:86
Vec_Z
vector< double > Vec_Z
Definition: CgemLineFit.cxx:85
f20
int f20(32)
Vec_m_phi
vector< double > Vec_m_phi
Definition: CgemLineFit.cxx:86
Mp
CgemMidDriftPlane * Mp
Definition: CgemLineFit.cxx:88
fa
int fa(63)
Vec_m_Q
vector< double > Vec_m_Q
Definition: CgemLineFit.cxx:86
Vec_m_Z
vector< double > Vec_m_Z
Definition: CgemLineFit.cxx:86
f00
int f00(4)
length
double length
Definition: CgemLineFit.cxx:106
GeoLayer0
CgemGeoLayer * GeoLayer0
Definition: CgemLineFit.cxx:90
Vec_sheetid
vector< double > Vec_sheetid
Definition: CgemLineFit.cxx:85
f10
int f10(16)
Vec_m_v
vector< double > Vec_m_v
Definition: CgemLineFit.cxx:86
l_outer
StraightLine * l_outer
Definition: CgemLineFit.cxx:93
pl_00
CgemGeoReadoutPlane * pl_00
Definition: CgemLineFit.cxx:95
_clst_0
int _clst_0(0)
tanl_set
double tanl_set
Definition: CgemLineFit.cxx:106
NC
int NC
Definition: CgemLineFit.cxx:108
GeoLayer1
CgemGeoLayer * GeoLayer1
Definition: CgemLineFit.cxx:91
Vec_layer
vector< double > Vec_layer
Definition: CgemLineFit.cxx:85
Vec_m_flag
vector< double > Vec_m_flag
Definition: CgemLineFit.cxx:86
Vec_v
vector< double > Vec_v
Definition: CgemLineFit.cxx:85
_clst_1
int _clst_1(0)
pl_20
CgemGeoReadoutPlane * pl_20
Definition: CgemLineFit.cxx:98
dz_set
double dz_set
Definition: CgemLineFit.cxx:106
Align_FLAG
bool Align_FLAG(false)
Vec_phi
vector< double > Vec_phi
Definition: CgemLineFit.cxx:85
sheet_flag
int sheet_flag(0)
Vec_layerid
vector< double > Vec_layerid
Definition: CgemLineFit.cxx:85
dr_set
double dr_set
Definition: CgemLineFit.cxx:106
pl_11
CgemGeoReadoutPlane * pl_11
Definition: CgemLineFit.cxx:97
pl_10
CgemGeoReadoutPlane * pl_10
Definition: CgemLineFit.cxx:96
pl_21
CgemGeoReadoutPlane * pl_21
Definition: CgemLineFit.cxx:99
Vec_m_sheetid
vector< double > Vec_m_sheetid
Definition: CgemLineFit.cxx:86
phi2
Double_t phi2
Definition: DataBase/tau_mode.c:7
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
phi1
Double_t phi1
Definition: DataBase/tau_mode.c:7
arg
double arg(const EvtComplex &c)
Definition: EvtComplex.hh:227
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
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
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: RecMdcHit.h:99
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: RecMdcTrack.h:102
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
Track
#define Track
Definition: TTrackBase.h:30
_Z
#define _Z
Definition: cfortran.h:871
CgemGeoAlign::getDx
double getDx(int layer_vir)
Definition: CgemGeoAlign.h:69
CgemGeoAlign::getRx
double getRx(int layer_vir)
Definition: CgemGeoAlign.h:72
CgemGeoAlign::getRy
double getRy(int layer_vir)
Definition: CgemGeoAlign.h:73
CgemGeoAlign::getRz
double getRz(int layer_vir)
Definition: CgemGeoAlign.h:74
CgemGeoAlign::getDz
double getDz(int layer_vir)
Definition: CgemGeoAlign.h:71
CgemGeoAlign::point_invTransform
HepPoint3D point_invTransform(int layer_vir, HepPoint3D pos)
Definition: CgemGeoAlign.cxx:322
CgemGeoAlign::getDy
double getDy(int layer_vir)
Definition: CgemGeoAlign.h:70
CgemGeoLayer::getMiddleROfGapD
double getMiddleROfGapD() const
Definition: CgemGeoLayer.h:105
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition: CgemGeomSvc.h:48
CgemGeomSvc::getLengthOfCgem
double getLengthOfCgem() const
Definition: CgemGeomSvc.h:42
CgemGeomSvc::getReadoutPlane
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition: CgemGeomSvc.h:140
CgemGeomSvc::getMidDriftPtr
CgemMidDriftPlane * getMidDriftPtr() const
Definition: CgemGeomSvc.h:72
CgemGeomSvc::getAlignPtr
CgemGeoAlign * getAlignPtr() const
Definition: CgemGeomSvc.h:69
CgemLineFit::Rec_Clusters_mTPC
void Rec_Clusters_mTPC()
Definition: CgemLineFit.cxx:2559
CgemLineFit::Loop_All
bool Loop_All()
Definition: CgemLineFit.cxx:925
CgemLineFit::GetNMaxQ
static vector< int > GetNMaxQ(vector< double > Q_11, vector< int > L_11, int Nmax)
Definition: CgemLineFit.cxx:2924
CgemLineFit::finalize
StatusCode finalize()
Definition: CgemLineFit.cxx:577
CgemLineFit::OrderClusters
void OrderClusters()
Definition: CgemLineFit.cxx:2061
CgemLineFit::Filter
void Filter(int layerid, StraightLine *l1)
Definition: CgemLineFit.cxx:2189
CgemLineFit::erase
void erase(int i)
Definition: CgemLineFit.cxx:2130
CgemLineFit::Swap
void Swap(int i, int j)
Definition: CgemLineFit.cxx:2084
CgemLineFit::IniPar
StraightLine * IniPar(double phi1, double z1, int i, double phi2, double z2, int j)
Definition: CgemLineFit.cxx:2171
CgemLineFit::Get4Par
vector< double > Get4Par(HepLorentzVector p4, Hep3Vector bp)
Definition: CgemLineFit.cxx:2400
CgemLineFit::Data_Max
bool Data_Max()
Definition: CgemLineFit.cxx:656
CgemLineFit::Min_dis_up
double Min_dis_up(int R_id, double x, double z, StraightLine *l1)
Definition: CgemLineFit.cxx:2853
CgemLineFit::execute
StatusCode execute()
Definition: CgemLineFit.cxx:399
CgemLineFit::CgemLineFit
CgemLineFit(const std::string &name, ISvcLocator *pSvcLocator)
Definition: CgemLineFit.cxx:111
CgemLineFit::GetVirLay
int GetVirLay(int geolay, double phi)
Definition: CgemLineFit.cxx:2917
CgemLineFit::Fit_Clusters
void Fit_Clusters(double par[])
Definition: CgemLineFit.cxx:2610
CgemLineFit::Erase_outer
bool Erase_outer()
Definition: CgemLineFit.cxx:2105
CgemLineFit::fcn2
static void fcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: CgemLineFit.cxx:1877
CgemLineFit::Layer_cross
bool Layer_cross(int R_id, StraightLine *l1)
Definition: CgemLineFit.cxx:2888
CgemLineFit::Trans
static vector< double > Trans(double par0, double par1, double par2, double par3)
Definition: CgemLineFit.cxx:2904
CgemLineFit::Min_dis_down
double Min_dis_down(int R_id, double x, double z, StraightLine *l1)
Definition: CgemLineFit.cxx:2870
CgemLineFit::fcn3
static void fcn3(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: CgemLineFit.cxx:1898
CgemLineFit::fcn
static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: CgemLineFit.cxx:1859
CgemLineFit::IncidentAngle
void IncidentAngle(StraightLine *cosmic_ray, HepPoint3D cluster, double theta, double ang[])
Definition: CgemLineFit.cxx:3039
CgemLineFit::dV
static double dV(int layer, double R, double dr, double phi0, double z0, double tglam, double x, double V)
Definition: CgemLineFit.cxx:1951
CgemLineFit::Loop_MaxQ
bool Loop_MaxQ()
Definition: CgemLineFit.cxx:1351
CgemLineFit::Min_dis
double Min_dis(int R_id, double x, double z, StraightLine *l1)
Definition: CgemLineFit.cxx:2830
CgemLineFit::MC_truth
vector< double > MC_truth()
Definition: CgemLineFit.cxx:2360
CgemLineFit::HitPos
bool HitPos(HepLorentzVector p4, Hep3Vector bp)
Definition: CgemLineFit.cxx:2425
CgemLineFit::Switch_CCmTPC
int Switch_CCmTPC
Definition: CgemLineFit.h:65
CgemLineFit::to_0_2pi
static void to_0_2pi(double &arg)
Definition: CgemLineFit.cxx:2814
CgemLineFit::Data_Closest
bool Data_Closest()
Definition: CgemLineFit.cxx:814
CgemLineFit::Fit
TMinuit * Fit(StraightLine *l_outer, int sheetflag, int typ)
Definition: CgemLineFit.cxx:1971
CgemLineFit::dx
static double dx(double layerid, double R, double dr, double phi0, double z0, double tanl, double x)
Definition: CgemLineFit.cxx:1924
CgemLineFit::Get_OtherIniPar
void Get_OtherIniPar(int clst_0, int clst_1, int clst_2)
Definition: CgemLineFit.cxx:2453
CgemLineFit::GetIntersection
void GetIntersection()
Definition: CgemLineFit.cxx:2637
CgemLineFit::FilterClusters
void FilterClusters()
Definition: CgemLineFit.cxx:2964
CgemLineFit::Get_MCInf
bool Get_MCInf()
Definition: CgemLineFit.cxx:2326
CgemLineFit::Min_diff
static double Min_diff(double arg1, double arg2)
Definition: CgemLineFit.cxx:2821
CgemLineFit::initialize
StatusCode initialize()
Definition: CgemLineFit.cxx:130
CgemLineFit::Store_Trk
void Store_Trk(TMinuit *fit, int trkid)
Definition: CgemLineFit.cxx:2987
CgemLineFit::Discard
void Discard(int layer)
Definition: CgemLineFit.cxx:2163
CgemLineFit::MinQ_Clus2D
double MinQ_Clus2D
Definition: CgemLineFit.h:110
CgemLineFit::RealPhi
double RealPhi(double SimuPhi)
Definition: CgemLineFit.cxx:3138
CgemLineFit::Get_Clusters
vector< double > Get_Clusters(vector< double >Vec_truth)
Definition: CgemLineFit.cxx:2759
CgemLineFit::Rec_Clusters
void Rec_Clusters()
Definition: CgemLineFit.cxx:2507
CgemLineFit::registerTrack
StatusCode registerTrack(RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
Definition: CgemLineFit.cxx:1806
CgemMidDriftPlane::getPointIdealGeom
bool getPointIdealGeom(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
Definition: CgemMidDriftPlane.cxx:24
CgemMidDriftPlane::getPointAligned
bool getPointAligned(int layer_vir, StraightLine pLine, HepPoint3D &posUp, HepPoint3D &posDown, double phiVUp[], double phiVDown[])
Definition: CgemMidDriftPlane.cxx:164
CgemMidDriftPlane::getR
double getR(int layer)
Definition: CgemMidDriftPlane.h:19
DstMdcTrack::setTrackId
void setTrackId(const int trackId)
Definition: DstMdcTrack.h:84
DstMdcTrack::setError
void setError(double err[15])
Definition: DstMdcTrack.cxx:118
DstMdcTrack::setHelix
void setHelix(double helix[5])
Definition: DstMdcTrack.cxx:104
DstMdcTrack::setChi2
void setChi2(const double chi)
Definition: DstMdcTrack.h:98
RecCgemCluster::getrecphi_mTPC
double getrecphi_mTPC(void) const
Definition: RecCgemCluster.h:65
RecCgemCluster::getenergydeposit
double getenergydeposit(void) const
Definition: RecCgemCluster.h:58
RecCgemCluster::getRecZ
double getRecZ(void) const
Definition: RecCgemCluster.h:61
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition: RecCgemCluster.h:53
RecCgemCluster::getrecv_mTPC
double getrecv_mTPC(void) const
Definition: RecCgemCluster.h:66
RecCgemCluster::getlayerid
int getlayerid(void) const
Definition: RecCgemCluster.h:55
RecCgemCluster::getflag
int getflag(void) const
Definition: RecCgemCluster.h:57
RecCgemCluster::getclusterflagb
int getclusterflagb(void) const
Definition: RecCgemCluster.h:71
RecCgemCluster::getRecZ_mTPC
double getRecZ_mTPC(void) const
Definition: RecCgemCluster.h:67
RecCgemCluster::getrecphi
double getrecphi(void) const
Definition: RecCgemCluster.h:59
RecCgemCluster::getrecv
double getrecv(void) const
Definition: RecCgemCluster.h:60
RecCgemCluster::setTrkId
void setTrkId(int trkid)
Definition: RecCgemCluster.h:36
RecCgemCluster::getsheetid
int getsheetid(void) const
Definition: RecCgemCluster.h:56
RecCgemCluster::getclusterflage
int getclusterflage(void) const
Definition: RecCgemCluster.h:72
RecMdcTrack::setVecClusters
void setVecClusters(ClusterRefVec vecclusters)
Definition: RecMdcTrack.cxx:86
RecMdcTrack::setNcluster
void setNcluster(int ncluster)
Definition: RecMdcTrack.h:83
StraightLine::dz
double dz(void) const
Definition: StraightLine.h:48
StraightLine::x
HepPoint3D x(double s=0.) const
returns position after moving s in downwoards
Definition: StraightLine.cxx:81
StraightLine::dr
double dr(void) const
returns an element of parameters.
Definition: StraightLine.h:45
StraightLine::tanl
double tanl(void) const
Definition: StraightLine.h:49
StraightLine::phi0
double phi0(void) const
Definition: StraightLine.h:46
H
IMPLICIT REAL *A H
Definition: myXsection.h:1
Event
Definition: Event.h:21
t
int t()
Definition: t.c:1