CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemLineFit.h
Go to the documentation of this file.
1#ifndef CGEMLINEFIT_H
2#define CGEMLINEFIT_H
3
4#include "GaudiKernel/Algorithm.h"
5#include "GaudiKernel/NTuple.h"
6#include "GaudiKernel/INTupleSvc.h"
7#include "MdcGeomSvc/IMdcGeomSvc.h"
8#include "MdcGeomSvc/MdcGeoWire.h"
9#include "MdcGeomSvc/MdcGeoLayer.h"
13#include <string>
14#include <vector>
15#include <map>
16#include "TH2D.h"
17#include "TF1.h"
18#include "TGraph.h"
19#include <fstream>
20
21#include "CLHEP/Alist/AList.h"
22#include "MdcGeom/MdcDetector.h"
23#include "BField/BField.h"
29//#include "CLHEP/Matrix/Vector.h"
30
31
32
35#include "HepPDT/ParticleDataTable.hh"
39#include "TrkBase/TrkRecoTrk.h"
43#include "TrackUtil/Helix.h"
44#include "TMinuit.h"
45class IBesTimerSvc;
46class BesTimer;
47
48
49class CgemLineFit:public Algorithm {
50 public:
51 CgemLineFit(const std::string& name, ISvcLocator* pSvcLocator);
53 StatusCode initialize();
54 StatusCode execute();
55 StatusCode finalize();
56
57
58 bool Data_Max();
59 bool Data_Closest();
60 bool Loop_All();
61 bool Loop_MaxQ();
62 bool ToyMC();
63 bool debug;
64 bool MC;
66 int TEST_N;
67 int Nmax;// the number of max cluster No. on eacy sheet
69 void OrderClusters ();
70 void OrderClusterSizeQ ();
71 void FilterClusters ();
72 bool Erase_outer();
73 void erase (int i);
74 int GetVirLay(int geolay, double phi);
75 void Discard (int layer);
76 StraightLine* IniPar(double phi1,double z1,int i,double phi2,double z2,int j);
77 void Filter (int layerid,StraightLine* l1);
78 static void fcn( Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
79 static void fcn2( Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
80 static void fcn3( Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
81 void Store_Trk(TMinuit* fit,int trkid);
82 TMinuit* Fit(StraightLine* l_outer,int sheetflag,int typ );
83
84 static double dx(double layerid, double R,double dr,double phi0,double z0,double tanl,double x );
85 static double dV(int layer,double R ,double dr,double phi0,double z0,double tglam,double x,double V);
86 static double CalSigma2(double layerid ,double flag,double x,double V,double z,StraightLine* cosmic_ray);
87 static double CalLineSAtY(StraightLine& l,double y);
88 bool Get_MCInf();
89 static void IncidentAngle(StraightLine* cosmic_ray,HepPoint3D cluster,double theta,double ang[] );
90 void GetIntersection(StraightLine* line);
91 vector<double>MC_truth();
92 vector<double>Get4Par(HepLorentzVector p4,Hep3Vector bp);
93 bool HitPos(HepLorentzVector p4,Hep3Vector bp);
94 void Get_OtherIniPar(int clst_0,int clst_1,int clst_2);
95 vector<double>Get_Clusters(vector<double>Vec_truth);
96 void Rec_Clusters();
97 void Rec_Clusters_mTPC();
98 void Rec_ClusterSize();
99 void Rec_ClusterQ();
100 void fireStripQ();
101 void fireStripT();
102 void fireStripID();
103 void Fit_Clusters(double par[]);
105 double RealPhi(double SimuPhi);
106
107 static vector<double> Trans(double par0,double par1,double par2,double par3);
108 static vector<int> GetNMaxQ(vector<double> Q_11,vector<int> L_11,int Nmax);
109 static double func(double layer,double R,double dr,double phi0,double z0,double tglam,double V,double x);
110
111 StatusCode registerTrack(RecMdcTrackCol*& recMdcTrackCol ,RecMdcHitCol*& recMdcHitCol);
112 void Erase(int _layerid);
113 static void to_0_2pi(double &arg);
114 static double Min_diff(double arg1,double arg2);
115 //static double Min_diff2(double arg1,double arg2);
116
117
121
122 double Min_dis(int R_id,double x,double z,StraightLine* l1 );
123 double Min_dis_up(int R_id,double x,double z,StraightLine* l1 );
124 double Min_dis_down(int R_id,double x,double z,StraightLine* l1 );
125 bool Layer_cross(int R_id,StraightLine* l1 );
126 void Swap(int i,int j);
127 private:
128
129 ClusterRefVec vecclusters;
130 double length;
131 double sigma;
132 double Chisq_cut;
133 double MAX_Clusters;
134 bool Align_Flag;
135 bool MC_Flag;
136 bool Run10_Flag;
137 string m_LUTfile;
138 ICgemCalibFunSvc* myCgemCalibSvc;
139 CgemLUTReader *lutreader;
140
141 bool ThreeClsts_Flag;
142 int Method;
143 RecMdcHitCol* m_hitList_tds;
144
145 float get_Time(CgemDigiCol::iterator iDigiCol);
146 float get_TimeWalk(CgemDigiCol::iterator iDigiCol);
147 float get_TimeRising(CgemDigiCol::iterator iDigiCol);
148 float get_TimeReference(CgemDigiCol::iterator iDigiCol);
149
150 NTuple::Tuple* m_tuple_track;
151 NTuple::Item<int> tr_run;
152 NTuple::Item<int> tr_event;
153 // track parameters
154 NTuple::Item<double> tr_drho;
155 NTuple::Item<double> tr_phi0;
156 NTuple::Item<double> tr_dz0;
157 NTuple::Item<double> tr_tanLam;
158 NTuple::Item<double> tr_chisq;
159 // whether is the track fitted successfully
160 NTuple::Item<bool> tr_Is_fitted;
161 // number of clusters in an event
162 NTuple::Item<int> tr_ncluster;
163
164 // cluster id
165 NTuple::Array<int> tr_id_cluster;
166 // cluster position
167 NTuple::Array<double> tr_x_cluster;
168 NTuple::Array<double> tr_y_cluster;
169 NTuple::Array<double> tr_z_cluster;
170 NTuple::Array<double> tr_phi_cluster;
171 // cluster deposit energy, layer id, sheet id and whether is it one selected by track on each sheet
172 NTuple::Array<double> tr_Q_cluster;
173 NTuple::Array<int> tr_layer_cluster;
174 NTuple::Array<int> tr_sheet_cluster;
175 NTuple::Array<int> tr_Is_selected_cluster;
176
177 NTuple::Tuple* m_tuple;
178 NTuple::Item<double> run;
179
180 NTuple::Item<double> R0;
181 NTuple::Item<double> R1;
182 NTuple::Item<double> R2;
183
184 NTuple::Item<double> status1;
185 NTuple::Item<double> status2;
186 NTuple::Item<double> status3;
187
188 NTuple::Item<double> clst_2;
189 NTuple::Item<double> clst_0;
190 NTuple::Item<double> clst_1;
191
192 NTuple::Item<double> clst_00;
193 NTuple::Item<double> clst_01;
194
195 NTuple::Item<double> clst_10;
196 NTuple::Item<double> clst_11;
197
198 NTuple::Item<double> clst_20;
199 NTuple::Item<double> clst_21;
200
201 NTuple::Item<double> event;
202 NTuple::Item<double> phi_o1;
203 NTuple::Item<double> phi_o2;
204 NTuple::Item<double> ini_1;
205 NTuple::Item<double> ini_0;
206 NTuple::Item<double> ini_2;
207 NTuple::Item<double> ini_3;
208
209
210 NTuple::Item<double> inim_1;
211 NTuple::Item<double> inim_0;
212 NTuple::Item<double> inim_2;
213 NTuple::Item<double> inim_3;
214 NTuple::Item<double> inii_1;
215 NTuple::Item<double> inii_0;
216 NTuple::Item<double> inii_2;
217 NTuple::Item<double> inii_3;
218 NTuple::Item<double> par0;
219 NTuple::Item<double> par1;
220 NTuple::Item<double> par2;
221 NTuple::Item<double> par3;
222 NTuple::Item<double> Err_par0;
223 NTuple::Item<double> Err_par1;
224 NTuple::Item<double> Err_par2;
225 NTuple::Item<double> Err_par3;
226 NTuple::Item<double> newpar0;
227 NTuple::Item<double> newpar1;
228 NTuple::Item<double> newpar2;
229 NTuple::Item<double> newpar3;
230 NTuple::Item<double> newErr_par0;
231 NTuple::Item<double> newErr_par1;
232 NTuple::Item<double> newErr_par2;
233 NTuple::Item<double> newErr_par3;
234 NTuple::Item<double> MC_par0;
235 NTuple::Item<double> MC_par1;
236 NTuple::Item<double> MC_par2;
237 NTuple::Item<double> MC_par3;
238 NTuple::Item<double> MC_px;
239 NTuple::Item<double> MC_py;
240 NTuple::Item<double> MC_pz;
241
242 NTuple::Item<double> x_0_up;
243 NTuple::Item<double> x_1_up;
244 NTuple::Item<double> x_2_up;
245 NTuple::Item<double> x_0_down;
246 NTuple::Item<double> x_1_down;
247 NTuple::Item<double> x_2_down;
248 NTuple::Item<double> v_0_up;
249 NTuple::Item<double> v_1_up;
250 NTuple::Item<double> v_2_up;
251 NTuple::Item<double> v_0_down;
252 NTuple::Item<double> v_1_down;
253 NTuple::Item<double> v_2_down;
254 NTuple::Item<double> z_1_up;
255 NTuple::Item<double> z_2_up;
256 NTuple::Item<double> z_0_up;
257 NTuple::Item<double> z_1_down;
258 NTuple::Item<double> z_2_down;
259 NTuple::Item<double> z_0_down;
260
261
262 NTuple::Item<double> x_0_up_m;
263 NTuple::Item<double> x_1_up_m;
264 NTuple::Item<double> x_2_up_m;
265 NTuple::Item<double> x_0_down_m;
266 NTuple::Item<double> x_1_down_m;
267 NTuple::Item<double> x_2_down_m;
268 NTuple::Item<double> v_0_up_m;
269 NTuple::Item<double> v_1_up_m;
270 NTuple::Item<double> v_2_up_m;
271 NTuple::Item<double> v_0_down_m;
272 NTuple::Item<double> v_1_down_m;
273 NTuple::Item<double> v_2_down_m;
274 NTuple::Item<double> z_1_up_m;
275 NTuple::Item<double> z_2_up_m;
276 NTuple::Item<double> z_0_up_m;
277 NTuple::Item<double> z_1_down_m;
278 NTuple::Item<double> z_2_down_m;
279 NTuple::Item<double> z_0_down_m;
280
281
282 NTuple::Item<double> x_0_up_f;
283 NTuple::Item<double> x_1_up_f;
284 NTuple::Item<double> x_2_up_f;
285 NTuple::Item<double> x_0_down_f;
286 NTuple::Item<double> x_1_down_f;
287 NTuple::Item<double> x_2_down_f;
288 NTuple::Item<double> v_0_up_f;
289 NTuple::Item<double> v_1_up_f;
290 NTuple::Item<double> v_2_up_f;
291 NTuple::Item<double> v_0_down_f;
292 NTuple::Item<double> v_1_down_f;
293 NTuple::Item<double> v_2_down_f;
294
295 NTuple::Item<double> x_0_up_mc;
296 NTuple::Item<double> x_1_up_mc;
297 NTuple::Item<double> x_2_up_mc;
298 NTuple::Item<double> x_0_down_mc;
299 NTuple::Item<double> x_1_down_mc;
300 NTuple::Item<double> x_2_down_mc;
301 NTuple::Item<double> v_0_up_mc;
302 NTuple::Item<double> v_1_up_mc;
303 NTuple::Item<double> v_2_up_mc;
304 NTuple::Item<double> v_0_down_mc;
305 NTuple::Item<double> v_1_down_mc;
306 NTuple::Item<double> v_2_down_mc;
307 NTuple::Item<double> x_0_up_mc2;
308 NTuple::Item<double> x_1_up_mc2;
309 NTuple::Item<double> x_2_up_mc2;
310 NTuple::Item<double> x_0_down_mc2;
311 NTuple::Item<double> x_1_down_mc2;
312 NTuple::Item<double> x_2_down_mc2;
313 NTuple::Item<double> v_0_up_mc2;
314 NTuple::Item<double> v_1_up_mc2;
315 NTuple::Item<double> v_2_up_mc2;
316 NTuple::Item<double> v_0_down_mc2;
317 NTuple::Item<double> v_1_down_mc2;
318 NTuple::Item<double> v_2_down_mc2;
319 NTuple::Item<double> z_0_up_mc2;
320 NTuple::Item<double> z_1_up_mc2;
321 NTuple::Item<double> z_2_up_mc2;
322 NTuple::Item<double> z_0_down_mc2;
323 NTuple::Item<double> z_1_down_mc2;
324 NTuple::Item<double> z_2_down_mc2;
325
326 NTuple::Item<int> x_0_up_size;
327 NTuple::Item<int> x_1_up_size;
328 NTuple::Item<int> x_2_up_size;
329 NTuple::Item<int> x_0_down_size;
330 NTuple::Item<int> x_1_down_size;
331 NTuple::Item<int> x_2_down_size;
332 NTuple::Item<int> v_0_up_size;
333 NTuple::Item<int> v_1_up_size;
334 NTuple::Item<int> v_2_up_size;
335 NTuple::Item<int> v_0_down_size;
336 NTuple::Item<int> v_1_down_size;
337 NTuple::Item<int> v_2_down_size;
338
339 NTuple::Array<int> x_0_up_stripID;
340 NTuple::Array<int> x_1_up_stripID;
341 NTuple::Array<int> x_2_up_stripID;
342 NTuple::Array<int> x_0_down_stripID;
343 NTuple::Array<int> x_1_down_stripID;
344 NTuple::Array<int> x_2_down_stripID;
345 NTuple::Array<int> v_0_up_stripID;
346 NTuple::Array<int> v_1_up_stripID;
347 NTuple::Array<int> v_2_up_stripID;
348 NTuple::Array<int> v_0_down_stripID;
349 NTuple::Array<int> v_1_down_stripID;
350 NTuple::Array<int> v_2_down_stripID;
351
352 NTuple::Item<double> x_0_up_Q;
353 NTuple::Item<double> x_1_up_Q;
354 NTuple::Item<double> x_2_up_Q;
355 NTuple::Item<double> x_0_down_Q;
356 NTuple::Item<double> x_1_down_Q;
357 NTuple::Item<double> x_2_down_Q;
358 NTuple::Item<double> v_0_up_Q;
359 NTuple::Item<double> v_1_up_Q;
360 NTuple::Item<double> v_2_up_Q;
361 NTuple::Item<double> v_0_down_Q;
362 NTuple::Item<double> v_1_down_Q;
363 NTuple::Item<double> v_2_down_Q;
364
365 NTuple::Array<double> x_0_up_stripQ;
366 NTuple::Array<double> x_1_up_stripQ;
367 NTuple::Array<double> x_2_up_stripQ;
368 NTuple::Array<double> x_0_down_stripQ;
369 NTuple::Array<double> x_1_down_stripQ;
370 NTuple::Array<double> x_2_down_stripQ;
371 NTuple::Array<double> v_0_up_stripQ;
372 NTuple::Array<double> v_1_up_stripQ;
373 NTuple::Array<double> v_2_up_stripQ;
374 NTuple::Array<double> v_0_down_stripQ;
375 NTuple::Array<double> v_1_down_stripQ;
376 NTuple::Array<double> v_2_down_stripQ;
377
378 NTuple::Array<double> x_0_up_stripT;
379 NTuple::Array<double> x_1_up_stripT;
380 NTuple::Array<double> x_2_up_stripT;
381 NTuple::Array<double> x_0_down_stripT;
382 NTuple::Array<double> x_1_down_stripT;
383 NTuple::Array<double> x_2_down_stripT;
384 NTuple::Array<double> v_0_up_stripT;
385 NTuple::Array<double> v_1_up_stripT;
386 NTuple::Array<double> v_2_up_stripT;
387 NTuple::Array<double> v_0_down_stripT;
388 NTuple::Array<double> v_1_down_stripT;
389 NTuple::Array<double> v_2_down_stripT;
390
391 NTuple::Array<double> x_0_up_stripTf;
392 NTuple::Array<double> x_1_up_stripTf;
393 NTuple::Array<double> x_2_up_stripTf;
394 NTuple::Array<double> x_0_down_stripTf;
395 NTuple::Array<double> x_1_down_stripTf;
396 NTuple::Array<double> x_2_down_stripTf;
397 NTuple::Array<double> v_0_up_stripTf;
398 NTuple::Array<double> v_1_up_stripTf;
399 NTuple::Array<double> v_2_up_stripTf;
400 NTuple::Array<double> v_0_down_stripTf;
401 NTuple::Array<double> v_1_down_stripTf;
402 NTuple::Array<double> v_2_down_stripTf;
403
404 NTuple::Item<double> inner_phi1;
405 NTuple::Item<double> inner_phi2;
406 NTuple::Item<double> inner_V1;
407 NTuple::Item<double> inner_V2;
408 NTuple::Item<double> inner_z1;
409 NTuple::Item<double> inner_z2;
410 NTuple::Item<double> chisq_1;
411 NTuple::Item<double> chisq_2;
412 NTuple::Item<double> chisq_3;
413
414 NTuple::Item<double> pos_x;
415 NTuple::Item<double> pos_y;
416 NTuple::Item<double> pos_z;
417
418 NTuple::Item<double> hit_x;
419 NTuple::Item<double> hit_y;
420 NTuple::Item<double> hit_z;
421
422 NTuple::Item<double> ang_1;
423 NTuple::Item<double> ang_1_x;
424 NTuple::Item<double> ang_1_v;
425
426 NTuple::Item<double> ang_2;
427 NTuple::Item<double> ang_2_x;
428 NTuple::Item<double> ang_2_v;
429
430 NTuple::Item<double> ang_3;
431 NTuple::Item<double> ang_3_x;
432 NTuple::Item<double> ang_3_v;
433
434 NTuple::Item<double> ang_4;
435 NTuple::Item<double> ang_4_x;
436 NTuple::Item<double> ang_4_v;
437
438 NTuple::Item<double> ang_5;
439 NTuple::Item<double> ang_5_x;
440 NTuple::Item<double> ang_5_v;
441
442 NTuple::Item<double> ang_6;
443 NTuple::Item<double> ang_6_x;
444 NTuple::Item<double> ang_6_v;
445
446 NTuple::Item<double> Angle;
447 NTuple::Item<double> ang_x;
448 NTuple::Item<double> ang_v;
449
450
451 NTuple::Item<double> inter_1_x;
452 NTuple::Item<double> inter_1_V;
453 NTuple::Item<int> inter_1_flag;
454 NTuple::Item<double> inter_1_z;
455
456 NTuple::Item<double> inter_2_x;
457 NTuple::Item<double> inter_2_V;
458 NTuple::Item<int> inter_2_flag;
459 NTuple::Item<double> inter_2_z;
460
461 NTuple::Item<double> inter_3_x;
462 NTuple::Item<double> inter_3_V;
463 NTuple::Item<int> inter_3_flag;
464 NTuple::Item<double> inter_3_z;
465
466 NTuple::Item<double> inter_4_x;
467 NTuple::Item<double> inter_4_V;
468 NTuple::Item<int> inter_4_flag;
469 NTuple::Item<double> inter_4_z;
470
471 NTuple::Item<double> inter_5_x;
472 NTuple::Item<double> inter_5_V;
473 NTuple::Item<int> inter_5_flag;
474 NTuple::Item<double> inter_5_z;
475
476 NTuple::Item<double> inter_6_x;
477 NTuple::Item<double> inter_6_V;
478 NTuple::Item<int> inter_6_flag;
479 NTuple::Item<double> inter_6_z;
480
481 NTuple::Item<double> inter_x;
482 NTuple::Item<double> inter_V;
483 NTuple::Item<int> inter_flag;
484 NTuple::Item<double> inter_z;
485
486 NTuple::Item<double> Rec_phi;
487 NTuple::Item<double> Rec_V;
488 NTuple::Item<double> Rec_Z;
489 NTuple::Item<double> Test_phi;
490 NTuple::Item<double> Test_V;
491 NTuple::Item<double> Test_Z;
492 void GetRes( StraightLine* l1,int i_inter);
493 //data member
494 int m_event;
495 int m_run;
496 RawDataProviderSvc* m_rawDataProviderSvc;
497 CgemGeomSvc* m_cgemGeomSvc;
498 int myNCgemLayers;
499 int myNSheets[4];
500 double myRXstrips[4];
501 double myRVstrips[4];
502 CgemGeoLayer* myCgemLayer[4];
503 // CgemGeomSvc* myCgemGeomSvc;
504 CgemCalibFunSvc* m_cgemCalibFunSvc;
505
506 // ReadoutPlane* readoutPlane;
507 RecMdcTrackCol* m_trackList_tds;
508 // CgemGeoReadoutPlane* readoutPlane;
509
510};
511#endif
StraightLine * l_outer
Double_t phi2
Double_t phi1
double arg(const EvtComplex &c)
ObjectVector< RecMdcHit > RecMdcHitCol
Definition RecMdcHit.h:99
SmartRefVector< RecCgemCluster > ClusterRefVec
Definition RecMdcTrack.h:28
ObjectVector< RecMdcTrack > RecMdcTrackCol
void Rec_Clusters_mTPC()
static vector< int > GetNMaxQ(vector< double > Q_11, vector< int > L_11, int Nmax)
StatusCode finalize()
void OrderClusters()
void Rec_ClusterQ()
void OrderClusterSizeQ()
void Filter(int layerid, StraightLine *l1)
void InAngle(StraightLine sl)
void erase(int i)
static double CalSigma2(double layerid, double flag, double x, double V, double z, StraightLine *cosmic_ray)
void Swap(int i, int j)
StraightLine * IniPar(double phi1, double z1, int i, double phi2, double z2, int j)
vector< double > Get4Par(HepLorentzVector p4, Hep3Vector bp)
double Min_dis_up(int R_id, double x, double z, StraightLine *l1)
StatusCode execute()
CgemLineFit(const std::string &name, ISvcLocator *pSvcLocator)
int GetVirLay(int geolay, double phi)
void Fit_Clusters(double par[])
static void fcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
void GetIntersection(StraightLine *line)
bool Layer_cross(int R_id, StraightLine *l1)
static vector< double > Trans(double par0, double par1, double par2, double par3)
double Min_dis_down(int R_id, double x, double z, StraightLine *l1)
static void fcn3(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static void IncidentAngle(StraightLine *cosmic_ray, HepPoint3D cluster, double theta, double ang[])
static double dV(int layer, double R, double dr, double phi0, double z0, double tglam, double x, double V)
double Min_dis(int R_id, double x, double z, StraightLine *l1)
vector< double > MC_truth()
bool HitPos(HepLorentzVector p4, Hep3Vector bp)
int Switch_CCmTPC
Definition CgemLineFit.h:68
static double func(double layer, double R, double dr, double phi0, double z0, double tglam, double V, double x)
static void to_0_2pi(double &arg)
bool Data_Closest()
TMinuit * Fit(StraightLine *l_outer, int sheetflag, int typ)
static double dx(double layerid, double R, double dr, double phi0, double z0, double tanl, double x)
void Get_OtherIniPar(int clst_0, int clst_1, int clst_2)
void FilterClusters()
static double Min_diff(double arg1, double arg2)
StatusCode initialize()
void Erase(int _layerid)
void Store_Trk(TMinuit *fit, int trkid)
static double CalLineSAtY(StraightLine &l, double y)
void Discard(int layer)
void Rec_ClusterSize()
double MinQ_Clus2D
double RealPhi(double SimuPhi)
vector< double > Get_Clusters(vector< double >Vec_truth)
void Rec_Clusters()
StatusCode registerTrack(RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)