CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemLineFit.cxx
Go to the documentation of this file.
1#include "CgemLineFit/CgemLineFit.h"
2#include "EventModel/EventHeader.h"
3#include "GaudiKernel/AlgFactory.h"
4#include "GaudiKernel/IDataProviderSvc.h"
5#include "GaudiKernel/INTupleSvc.h"
6#include "GaudiKernel/IPartPropSvc.h"
7#include "GaudiKernel/IService.h"
8#include "GaudiKernel/ISvcLocator.h"
9#include "GaudiKernel/MsgStream.h"
10#include "GaudiKernel/NTuple.h"
11#include "GaudiKernel/PropertyMgr.h"
12#include "GaudiKernel/SmartDataPtr.h"
13#include "Identifier/Identifier.h"
14#include "Identifier/MdcID.h"
15#include "MdcRawEvent/MdcDigi.h"
16#include "TStopwatch.h"
17#include <iostream>
18#include <math.h>
19#include <vector>
20
21#include "EvTimeEvent/RecEsTime.h"
22#include "MdcData/MdcHit.h"
23#include "MdcGeom/EntranceAngle.h"
24#include "RawEvent/RawDataUtil.h"
25
26#include "McTruth/CgemMcHit.h"
27#include "McTruth/DecayMode.h"
28#include "McTruth/McParticle.h"
29#include "MdcRecoUtil/Pdt.h"
30#include "TrackUtil/Helix.h"
31
32#include "MdcData/MdcRecoHitOnTrack.h"
33#include "MdcPrintSvc/IMdcPrintSvc.h"
34#include "MdcxReco/Mdcxprobab.h"
35#include "TrkBase/TrkExchangePar.h"
36#include "TrkBase/TrkFit.h"
37#include "TrkBase/TrkHitList.h"
38#include "TrkFitter/TrkCircleMaker.h"
39#include "TrkFitter/TrkHelixFitter.h"
40#include "TrkFitter/TrkHelixMaker.h"
41#include "TrkFitter/TrkLineMaker.h"
42
43#include <math.h>
44
45#include "CgemCalibFunSvc/CgemCalibFunSvc.h"
46#include "CgemGeomSvc/CgemGeomSvc.h"
47#include "CgemGeomSvc/ICgemGeomSvc.h"
48#include "EvTimeEvent/RecEsTime.h"
49#include "GaudiKernel/IPartPropSvc.h"
50#include "McTruth/MdcMcHit.h"
51#include "MdcData/MdcHit.h"
52#include "MdcGeom/EntranceAngle.h"
53#include "MdcPrintSvc/IMdcPrintSvc.h"
54#include "MdcRecoUtil/Pdt.h"
55#include "MdcTrkRecon/MdcMap.h"
56#include "RawEvent/RawDataUtil.h"
57#include "TrackUtil/Helix.h"
58#include "TrkBase/TrkExchangePar.h"
59#include "TrkBase/TrkFit.h"
60#include "TrkBase/TrkHitList.h"
61#include "TrkFitter/TrkCircleMaker.h"
62#include "TrkFitter/TrkHelixFitter.h"
63#include "TrkFitter/TrkHelixMaker.h"
64#include "TrkFitter/TrkLineMaker.h"
65
66#include "TH2D.h"
67#include "TTree.h"
68
69#include "BesTimerSvc/BesTimer.h"
70#include "BesTimerSvc/BesTimerSvc.h"
71#include "BesTimerSvc/IBesTimerSvc.h"
72#include "CgemRawEvent/CgemDigi.h"
73#include "Identifier/CgemID.h"
74#include "TArrayI.h"
75#include "TF1.h"
76#include "TGraph.h"
77#include "TGraphErrors.h"
78#include "TMath.h"
79#include "TMinuit.h"
80#include "TRandom.h"
81
82#include "CLHEP/Vector/LorentzVector.h" //
83#include "CLHEP/Vector/ThreeVector.h"
84#include "CLHEP/Vector/TwoVector.h"
85using CLHEP::Hep2Vector;
86using CLHEP::Hep3Vector;
87using CLHEP::HepLorentzVector;
88using namespace std;
89using namespace Event;
90
91vector<double> Vec_layer, Vec_phi, Vec_Z, Vec_layerid, Vec_Virlayerid, Vec_v, Vec_flag, Vec_Q, Vec_sheetid, Vec_xCluSize,
92 Vec_vCluSize, Vec_XQ, Vec_VQ; //
93vector<double> Vec_x_stripsQ, Vec_v_stripsQ, Vec_x_stripsT, Vec_v_stripsT, Vec_x_stripsTf, Vec_v_stripsTf;
94vector<vector<double> > Vec_XstripQ, Vec_VstripQ, Vec_XstripT, Vec_VstripT, Vec_XstripTf, Vec_VstripTf;
95vector<int> Vec_x_stripsid, Vec_v_stripsid;
96vector<vector<int> > Vec_XstripID, Vec_VstripID;
97vector<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; //
98
99CgemMidDriftPlane *Mp;
100CgemGeoAlign *Al;
101CgemGeoLayer *GeoLayer0;
102CgemGeoLayer *GeoLayer1;
103CgemGeoLayer *GeoLayer2;
104StraightLine *l_outer;
105
106CgemGeoReadoutPlane *pl_00;
107CgemGeoReadoutPlane *pl_10;
108CgemGeoReadoutPlane *pl_11;
109CgemGeoReadoutPlane *pl_20;
110CgemGeoReadoutPlane *pl_21;
111
112// int is_cluster_exsited[6] = {0,0,0,0,0,0};
113// int nEvent_with_6cluster = 0;
114// int nEvent_with_4cluster = 0;
115
116bool Align_FLAG(false);
117double sigma2(0); // resoulution of hardware
118// flag for each sheet ?
119int f21(1), f11(2), f00(4), f01(8), f10(16), f20(32), fa(63), sheet_flag(0);
120
121double length, dr_set, phi0_set, dz_set, tanl_set;
122double R_layer[3];
123double AngleOfStereo[3];
124int NC;
125int _clst_0(0), _clst_1(0), _clst_2(0);
126
127CgemLineFit::CgemLineFit(const std::string &name, ISvcLocator *pSvcLocator) : Algorithm(name, pSvcLocator) {
128 declareProperty("sigma", sigma = 0.13); // unit:mm
129 declareProperty("Run10_Flag", Run10_Flag = false); // true:4 clusters, false: 3 clusters
130 declareProperty("Align_Flag", Align_Flag = false); // Alignment
131 declareProperty("Switch_CCmTPC", Switch_CCmTPC = 0); // 0:Center charge 1:mTPC 2: merged position
132 declareProperty("Method", Method = 0); // 0: ToyMC 1: max_charge 2:closest to intersection 3: Loop_all 4: Loop_MaxQ
133 declareProperty("MAX_COMB", MAX_COMB = 50); // The maximum number of combinations for method 3
134 declareProperty("MAX_Clusters", MAX_Clusters = 500); // The maximum number of combinations for method 3
135 declareProperty("MaxClu_Sheet", Nmax = 3); // Max number of clusters on each sheet for method 4
136 declareProperty("Chisq_cut", Chisq_cut = 2000); // Max chisq for the selected set of clusters
137 declareProperty("TEST_N", TEST_N = 0); // set the sheet you want to study the efficiency / resolution. 1
138 // to 6 from top to bottom. 0 use all the layers
139 declareProperty("Debug", debug = 0); // Switch of debug
140 declareProperty("MinQ_Clus2D", MinQ_Clus2D = 0); // Add charge cut on 2D cluster to reduce the noise
141 declareProperty("MC", MC = 0); // remove the hits on 3rd layer for MC
142 declareProperty("LUTfile", m_LUTfile = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_17.root");
143}
144CgemLineFit::~CgemLineFit() { delete lutreader; }
145
146// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
147StatusCode CgemLineFit::initialize() {
148 MsgStream log(msgSvc(), name());
149 sigma2 = sigma * sigma;
150 Align_FLAG = Align_Flag;
151 log << MSG::INFO << "in initialize()" << endreq;
152 StatusCode sc;
153 NTuplePtr cosmic(ntupleSvc(), "LineFit/cosmic");
154 if (cosmic) m_tuple = cosmic;
155 else {
156 m_tuple = ntupleSvc()->book("LineFit/cosmic", CLID_ColumnWiseTuple, "N-Tuple example");
157 if (m_tuple) {
158
159 sc = m_tuple->addItem("run", run);
160 sc = m_tuple->addItem("event", event);
161 sc = m_tuple->addItem("R0", R0);
162 sc = m_tuple->addItem("R1", R1);
163 sc = m_tuple->addItem("R2", R2);
164 sc = m_tuple->addItem("clst_0", clst_0);
165 sc = m_tuple->addItem("clst_1", clst_1);
166 sc = m_tuple->addItem("clst_2", clst_2);
167 sc = m_tuple->addItem("clst_01", clst_01);
168 sc = m_tuple->addItem("clst_00", clst_00);
169 sc = m_tuple->addItem("clst_11", clst_11);
170 sc = m_tuple->addItem("clst_10", clst_10);
171 sc = m_tuple->addItem("clst_21", clst_21);
172 sc = m_tuple->addItem("clst_20", clst_20);
173 sc = m_tuple->addItem("status1", status1);
174 sc = m_tuple->addItem("status2", status2);
175 sc = m_tuple->addItem("status3", status3);
176 sc = m_tuple->addItem("ini_0", ini_0);
177 sc = m_tuple->addItem("ini_1", ini_1);
178 sc = m_tuple->addItem("ini_2", ini_2);
179 sc = m_tuple->addItem("ini_3", ini_3);
180 sc = m_tuple->addItem("inim_0", inim_0);
181 sc = m_tuple->addItem("inim_1", inim_1);
182 sc = m_tuple->addItem("inim_2", inim_2);
183 sc = m_tuple->addItem("inim_3", inim_3);
184 sc = m_tuple->addItem("inii_0", inii_0);
185 sc = m_tuple->addItem("inii_1", inii_1);
186 sc = m_tuple->addItem("inii_2", inii_2);
187 sc = m_tuple->addItem("inii_3", inii_3);
188 sc = m_tuple->addItem("par0", par0);
189 sc = m_tuple->addItem("par1", par1);
190 sc = m_tuple->addItem("par2", par2);
191 sc = m_tuple->addItem("par3", par3);
192 sc = m_tuple->addItem("MC_par0", MC_par0);
193 sc = m_tuple->addItem("MC_par1", MC_par1);
194 sc = m_tuple->addItem("MC_par2", MC_par2);
195 sc = m_tuple->addItem("MC_par3", MC_par3);
196 sc = m_tuple->addItem("MC_px", MC_px);
197 sc = m_tuple->addItem("MC_py", MC_py);
198 sc = m_tuple->addItem("MC_pz", MC_pz);
199 sc = m_tuple->addItem("Err_par0", Err_par0);
200 sc = m_tuple->addItem("Err_par1", Err_par1);
201 sc = m_tuple->addItem("Err_par2", Err_par2);
202 sc = m_tuple->addItem("Err_par3", Err_par3);
203
204 sc = m_tuple->addItem("x_2_up", x_2_up);
205 sc = m_tuple->addItem("x_1_up", x_1_up);
206 sc = m_tuple->addItem("x_0_up", x_0_up);
207 sc = m_tuple->addItem("v_2_up", v_2_up);
208 sc = m_tuple->addItem("v_1_up", v_1_up);
209 sc = m_tuple->addItem("v_0_up", v_0_up);
210 sc = m_tuple->addItem("z_2_up", z_2_up);
211 sc = m_tuple->addItem("z_1_up", z_1_up);
212 sc = m_tuple->addItem("z_0_up", z_0_up);
213 sc = m_tuple->addItem("x_2_down", x_2_down);
214 sc = m_tuple->addItem("x_1_down", x_1_down);
215 sc = m_tuple->addItem("x_0_down", x_0_down);
216 sc = m_tuple->addItem("v_2_down", v_2_down);
217 sc = m_tuple->addItem("v_1_down", v_1_down);
218 sc = m_tuple->addItem("v_0_down", v_0_down);
219 sc = m_tuple->addItem("z_2_down", z_2_down);
220 sc = m_tuple->addItem("z_1_down", z_1_down);
221 sc = m_tuple->addItem("z_0_down", z_0_down);
222
223 sc = m_tuple->addItem("x_2_up_m", x_2_up_m);
224 sc = m_tuple->addItem("x_1_up_m", x_1_up_m);
225 sc = m_tuple->addItem("x_0_up_m", x_0_up_m);
226 sc = m_tuple->addItem("v_2_up_m", v_2_up_m);
227 sc = m_tuple->addItem("v_1_up_m", v_1_up_m);
228 sc = m_tuple->addItem("v_0_up_m", v_0_up_m);
229 sc = m_tuple->addItem("z_2_up_m", z_2_up_m);
230 sc = m_tuple->addItem("z_1_up_m", z_1_up_m);
231 sc = m_tuple->addItem("z_0_up_m", z_0_up_m);
232
233 sc = m_tuple->addItem("x_2_down_m", x_2_down_m);
234 sc = m_tuple->addItem("x_1_down_m", x_1_down_m);
235 sc = m_tuple->addItem("x_0_down_m", x_0_down_m);
236 sc = m_tuple->addItem("v_2_down_m", v_2_down_m);
237 sc = m_tuple->addItem("v_1_down_m", v_1_down_m);
238 sc = m_tuple->addItem("v_0_down_m", v_0_down_m);
239 sc = m_tuple->addItem("z_2_down_m", z_2_down_m);
240 sc = m_tuple->addItem("z_1_down_m", z_1_down_m);
241 sc = m_tuple->addItem("z_0_down_m", z_0_down_m);
242
243 sc = m_tuple->addItem("x_2_up_f", x_2_up_f);
244 sc = m_tuple->addItem("x_1_up_f", x_1_up_f);
245 sc = m_tuple->addItem("x_0_up_f", x_0_up_f);
246 sc = m_tuple->addItem("v_2_up_f", v_2_up_f);
247 sc = m_tuple->addItem("v_1_up_f", v_1_up_f);
248 sc = m_tuple->addItem("v_0_up_f", v_0_up_f);
249
250 sc = m_tuple->addItem("x_2_down_f", x_2_down_f);
251 sc = m_tuple->addItem("x_1_down_f", x_1_down_f);
252 sc = m_tuple->addItem("x_0_down_f", x_0_down_f);
253 sc = m_tuple->addItem("v_2_down_f", v_2_down_f);
254 sc = m_tuple->addItem("v_1_down_f", v_1_down_f);
255 sc = m_tuple->addItem("v_0_down_f", v_0_down_f);
256
257 sc = m_tuple->addItem("x_2_up_pre", x_2_up_mc);
258 sc = m_tuple->addItem("x_1_up_pre", x_1_up_mc);
259 sc = m_tuple->addItem("x_0_up_pre", x_0_up_mc);
260 sc = m_tuple->addItem("v_2_up_pre", v_2_up_mc);
261 sc = m_tuple->addItem("v_1_up_pre", v_1_up_mc);
262 sc = m_tuple->addItem("v_0_up_pre", v_0_up_mc);
263
264 sc = m_tuple->addItem("x_2_down_pre", x_2_down_mc);
265 sc = m_tuple->addItem("x_1_down_pre", x_1_down_mc);
266 sc = m_tuple->addItem("x_0_down_pre", x_0_down_mc);
267 sc = m_tuple->addItem("v_2_down_pre", v_2_down_mc);
268 sc = m_tuple->addItem("v_1_down_pre", v_1_down_mc);
269 sc = m_tuple->addItem("v_0_down_pre", v_0_down_mc);
270
271 sc = m_tuple->addItem("x_2_up_hit", x_2_up_mc2);
272 sc = m_tuple->addItem("x_1_up_hit", x_1_up_mc2);
273 sc = m_tuple->addItem("x_0_up_hit", x_0_up_mc2);
274 sc = m_tuple->addItem("v_2_up_hit", v_2_up_mc2);
275 sc = m_tuple->addItem("v_1_up_hit", v_1_up_mc2);
276 sc = m_tuple->addItem("v_0_up_hit", v_0_up_mc2);
277 sc = m_tuple->addItem("z_2_up_hit", z_2_up_mc2);
278 sc = m_tuple->addItem("z_1_up_hit", z_1_up_mc2);
279 sc = m_tuple->addItem("z_0_up_hit", z_0_up_mc2);
280
281 sc = m_tuple->addItem("x_2_down_hit", x_2_down_mc2);
282 sc = m_tuple->addItem("x_1_down_hit", x_1_down_mc2);
283 sc = m_tuple->addItem("x_0_down_hit", x_0_down_mc2);
284 sc = m_tuple->addItem("v_2_down_hit", v_2_down_mc2);
285 sc = m_tuple->addItem("v_1_down_hit", v_1_down_mc2);
286 sc = m_tuple->addItem("v_0_down_hit", v_0_down_mc2);
287 sc = m_tuple->addItem("z_2_down_hit", z_2_down_mc2);
288 sc = m_tuple->addItem("z_1_down_hit", z_1_down_mc2);
289 sc = m_tuple->addItem("z_0_down_hit", z_0_down_mc2);
290
291 sc = m_tuple->addItem("x_2_up_size", x_2_up_size, 0, 100);
292 sc = m_tuple->addItem("x_1_up_size", x_1_up_size, 0, 100);
293 sc = m_tuple->addItem("x_0_up_size", x_0_up_size, 0, 100);
294 sc = m_tuple->addItem("v_2_up_size", v_2_up_size, 0, 100);
295 sc = m_tuple->addItem("v_1_up_size", v_1_up_size, 0, 100);
296 sc = m_tuple->addItem("v_0_up_size", v_0_up_size, 0, 100);
297
298 sc = m_tuple->addItem("x_2_down_size", x_2_down_size, 0, 100);
299 sc = m_tuple->addItem("x_1_down_size", x_1_down_size, 0, 100);
300 sc = m_tuple->addItem("x_0_down_size", x_0_down_size, 0, 100);
301 sc = m_tuple->addItem("v_2_down_size", v_2_down_size, 0, 100);
302 sc = m_tuple->addItem("v_1_down_size", v_1_down_size, 0, 100);
303 sc = m_tuple->addItem("v_0_down_size", v_0_down_size, 0, 100);
304
305 sc = m_tuple->addItem("x_2_up_Q", x_2_up_Q);
306 sc = m_tuple->addItem("x_1_up_Q", x_1_up_Q);
307 sc = m_tuple->addItem("x_0_up_Q", x_0_up_Q);
308 sc = m_tuple->addItem("v_2_up_Q", v_2_up_Q);
309 sc = m_tuple->addItem("v_1_up_Q", v_1_up_Q);
310 sc = m_tuple->addItem("v_0_up_Q", v_0_up_Q);
311
312 sc = m_tuple->addItem("x_2_down_Q", x_2_down_Q);
313 sc = m_tuple->addItem("x_1_down_Q", x_1_down_Q);
314 sc = m_tuple->addItem("x_0_down_Q", x_0_down_Q);
315 sc = m_tuple->addItem("v_2_down_Q", v_2_down_Q);
316 sc = m_tuple->addItem("v_1_down_Q", v_1_down_Q);
317 sc = m_tuple->addItem("v_0_down_Q", v_0_down_Q);
318
319 sc = m_tuple->addIndexedItem("x_2_up_stripQ", x_2_up_size, x_2_up_stripQ);
320 sc = m_tuple->addIndexedItem("x_1_up_stripQ", x_1_up_size, x_1_up_stripQ);
321 sc = m_tuple->addIndexedItem("x_0_up_stripQ", x_0_up_size, x_0_up_stripQ);
322 sc = m_tuple->addIndexedItem("v_2_up_stripQ", v_2_up_size, v_2_up_stripQ);
323 sc = m_tuple->addIndexedItem("v_1_up_stripQ", v_1_up_size, v_1_up_stripQ);
324 sc = m_tuple->addIndexedItem("v_0_up_stripQ", v_0_up_size, v_0_up_stripQ);
325
326 sc = m_tuple->addIndexedItem("x_2_down_stripQ", x_2_down_size, x_2_down_stripQ);
327 sc = m_tuple->addIndexedItem("x_1_down_stripQ", x_1_down_size, x_1_down_stripQ);
328 sc = m_tuple->addIndexedItem("x_0_down_stripQ", x_0_down_size, x_0_down_stripQ);
329 sc = m_tuple->addIndexedItem("v_2_down_stripQ", v_2_down_size, v_2_down_stripQ);
330 sc = m_tuple->addIndexedItem("v_1_down_stripQ", v_1_down_size, v_1_down_stripQ);
331 sc = m_tuple->addIndexedItem("v_0_down_stripQ", v_0_down_size, v_0_down_stripQ);
332
333 sc = m_tuple->addIndexedItem("x_2_up_stripT", x_2_up_size, x_2_up_stripT);
334 sc = m_tuple->addIndexedItem("x_1_up_stripT", x_1_up_size, x_1_up_stripT);
335 sc = m_tuple->addIndexedItem("x_0_up_stripT", x_0_up_size, x_0_up_stripT);
336 sc = m_tuple->addIndexedItem("v_2_up_stripT", v_2_up_size, v_2_up_stripT);
337 sc = m_tuple->addIndexedItem("v_1_up_stripT", v_1_up_size, v_1_up_stripT);
338 sc = m_tuple->addIndexedItem("v_0_up_stripT", v_0_up_size, v_0_up_stripT);
339
340 sc = m_tuple->addIndexedItem("x_2_down_stripT", x_2_down_size, x_2_down_stripT);
341 sc = m_tuple->addIndexedItem("x_1_down_stripT", x_1_down_size, x_1_down_stripT);
342 sc = m_tuple->addIndexedItem("x_0_down_stripT", x_0_down_size, x_0_down_stripT);
343 sc = m_tuple->addIndexedItem("v_2_down_stripT", v_2_down_size, v_2_down_stripT);
344 sc = m_tuple->addIndexedItem("v_1_down_stripT", v_1_down_size, v_1_down_stripT);
345 sc = m_tuple->addIndexedItem("v_0_down_stripT", v_0_down_size, v_0_down_stripT);
346
347 sc = m_tuple->addIndexedItem("x_2_up_stripTf", x_2_up_size, x_2_up_stripTf);
348 sc = m_tuple->addIndexedItem("x_1_up_stripTf", x_1_up_size, x_1_up_stripTf);
349 sc = m_tuple->addIndexedItem("x_0_up_stripTf", x_0_up_size, x_0_up_stripTf);
350 sc = m_tuple->addIndexedItem("v_2_up_stripTf", v_2_up_size, v_2_up_stripTf);
351 sc = m_tuple->addIndexedItem("v_1_up_stripTf", v_1_up_size, v_1_up_stripTf);
352 sc = m_tuple->addIndexedItem("v_0_up_stripTf", v_0_up_size, v_0_up_stripTf);
353
354 sc = m_tuple->addIndexedItem("x_2_down_stripTf", x_2_down_size, x_2_down_stripTf);
355 sc = m_tuple->addIndexedItem("x_1_down_stripTf", x_1_down_size, x_1_down_stripTf);
356 sc = m_tuple->addIndexedItem("x_0_down_stripTf", x_0_down_size, x_0_down_stripTf);
357 sc = m_tuple->addIndexedItem("v_2_down_stripTf", v_2_down_size, v_2_down_stripTf);
358 sc = m_tuple->addIndexedItem("v_1_down_stripTf", v_1_down_size, v_1_down_stripTf);
359 sc = m_tuple->addIndexedItem("v_0_down_stripTf", v_0_down_size, v_0_down_stripTf);
360
361 sc = m_tuple->addIndexedItem("x_2_up_stripID", x_2_up_size, x_2_up_stripID);
362 sc = m_tuple->addIndexedItem("x_1_up_stripID", x_1_up_size, x_1_up_stripID);
363 sc = m_tuple->addIndexedItem("x_0_up_stripID", x_0_up_size, x_0_up_stripID);
364 sc = m_tuple->addIndexedItem("v_2_up_stripID", v_2_up_size, v_2_up_stripID);
365 sc = m_tuple->addIndexedItem("v_1_up_stripID", v_1_up_size, v_1_up_stripID);
366 sc = m_tuple->addIndexedItem("v_0_up_stripID", v_0_up_size, v_0_up_stripID);
367
368 sc = m_tuple->addIndexedItem("x_2_down_stripID", x_2_down_size, x_2_down_stripID);
369 sc = m_tuple->addIndexedItem("x_1_down_stripID", x_1_down_size, x_1_down_stripID);
370 sc = m_tuple->addIndexedItem("x_0_down_stripID", x_0_down_size, x_0_down_stripID);
371 sc = m_tuple->addIndexedItem("v_2_down_stripID", v_2_down_size, v_2_down_stripID);
372 sc = m_tuple->addIndexedItem("v_1_down_stripID", v_1_down_size, v_1_down_stripID);
373 sc = m_tuple->addIndexedItem("v_0_down_stripID", v_0_down_size, v_0_down_stripID);
374
375 sc = m_tuple->addItem("inter_1_x", inter_1_x);
376 sc = m_tuple->addItem("inter_2_x", inter_2_x);
377 sc = m_tuple->addItem("inter_3_x", inter_3_x);
378 sc = m_tuple->addItem("inter_4_x", inter_4_x);
379 sc = m_tuple->addItem("inter_5_x", inter_5_x);
380 sc = m_tuple->addItem("inter_6_x", inter_6_x);
381 sc = m_tuple->addItem("inter_x", inter_x);
382
383 sc = m_tuple->addItem("inter_1_V", inter_1_V);
384 sc = m_tuple->addItem("inter_2_V", inter_2_V);
385 sc = m_tuple->addItem("inter_3_V", inter_3_V);
386 sc = m_tuple->addItem("inter_4_V", inter_4_V);
387 sc = m_tuple->addItem("inter_5_V", inter_5_V);
388 sc = m_tuple->addItem("inter_6_V", inter_6_V);
389 sc = m_tuple->addItem("inter_V", inter_V);
390
391 sc = m_tuple->addItem("inter_1_flag", inter_1_flag);
392 sc = m_tuple->addItem("inter_2_flag", inter_2_flag);
393 sc = m_tuple->addItem("inter_3_flag", inter_3_flag);
394 sc = m_tuple->addItem("inter_4_flag", inter_4_flag);
395 sc = m_tuple->addItem("inter_5_flag", inter_5_flag);
396 sc = m_tuple->addItem("inter_6_flag", inter_6_flag);
397 sc = m_tuple->addItem("inter_flag", inter_flag);
398
399 sc = m_tuple->addItem("inter_1_z", inter_1_z);
400 sc = m_tuple->addItem("inter_2_z", inter_2_z);
401 sc = m_tuple->addItem("inter_3_z", inter_3_z);
402 sc = m_tuple->addItem("inter_4_z", inter_4_z);
403 sc = m_tuple->addItem("inter_5_z", inter_5_z);
404 sc = m_tuple->addItem("inter_6_z", inter_6_z);
405 sc = m_tuple->addItem("inter_z", inter_z);
406
407 sc = m_tuple->addItem("chisq_1", chisq_1);
408 sc = m_tuple->addItem("chisq_2", chisq_2);
409 sc = m_tuple->addItem("chisq_3", chisq_3);
410 sc = m_tuple->addItem("pos_x", pos_x);
411 sc = m_tuple->addItem("pos_y", pos_y);
412 sc = m_tuple->addItem("pos_z", pos_z);
413
414 sc = m_tuple->addItem("hit_x", hit_x);
415 sc = m_tuple->addItem("hit_y", hit_y);
416 sc = m_tuple->addItem("hit_z", hit_z);
417
418 sc = m_tuple->addItem("ang_1", ang_1);
419 sc = m_tuple->addItem("x_2_up_ang", ang_1_x);
420 sc = m_tuple->addItem("v_2_up_ang", ang_1_v);
421
422 sc = m_tuple->addItem("ang_2", ang_2);
423 sc = m_tuple->addItem("x_1_up_ang", ang_2_x);
424 sc = m_tuple->addItem("v_1_up_ang", ang_2_v);
425
426 sc = m_tuple->addItem("ang_3", ang_3);
427 sc = m_tuple->addItem("x_0_up_ang", ang_3_x);
428 sc = m_tuple->addItem("v_0_up_ang", ang_3_v);
429
430 sc = m_tuple->addItem("ang_4", ang_4);
431 sc = m_tuple->addItem("x_0_down_ang", ang_4_x);
432 sc = m_tuple->addItem("v_0_down_ang", ang_4_v);
433
434 sc = m_tuple->addItem("ang_5", ang_5);
435 sc = m_tuple->addItem("x_1_down_ang", ang_5_x);
436 sc = m_tuple->addItem("v_1_down_ang", ang_5_v);
437
438 sc = m_tuple->addItem("ang_6", ang_6);
439 sc = m_tuple->addItem("x_2_down_ang", ang_6_x);
440 sc = m_tuple->addItem("v_2_down_ang", ang_6_v);
441
442 sc = m_tuple->addItem("Angle", Angle);
443 sc = m_tuple->addItem("ang_x", ang_x);
444 sc = m_tuple->addItem("ang_v", ang_v);
445
446 sc = m_tuple->addItem("Rec_phi", Rec_phi);
447 sc = m_tuple->addItem("Rec_V", Rec_V);
448 sc = m_tuple->addItem("Rec_Z", Rec_Z);
449 sc = m_tuple->addItem("Test_phi", Test_phi);
450 sc = m_tuple->addItem("Test_V", Test_V);
451 sc = m_tuple->addItem("Test_Z", Test_Z);
452 } else {
453 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple) << endmsg;
454 return StatusCode::FAILURE;
455 }
456 }
457
458 NTuplePtr Track(ntupleSvc(), "LineFit/Track");
459 if (Track) m_tuple_track = Track;
460 else {
461 m_tuple_track = ntupleSvc()->book("LineFit/Track", CLID_ColumnWiseTuple, "N-Tuple example");
462 if (m_tuple_track) {
463 sc = m_tuple_track->addItem("run", tr_run);
464 sc = m_tuple_track->addItem("event", tr_event);
465 sc = m_tuple_track->addItem("drho", tr_drho);
466 sc = m_tuple_track->addItem("phi0", tr_phi0);
467 sc = m_tuple_track->addItem("dz0", tr_dz0);
468 sc = m_tuple_track->addItem("chisq", tr_chisq);
469 sc = m_tuple_track->addItem("Is_fitted", tr_Is_fitted);
470 sc = m_tuple_track->addItem("tanLam", tr_tanLam);
471 sc = m_tuple_track->addItem("ncluster", tr_ncluster, 0, 5000);
472 sc = m_tuple_track->addIndexedItem("id_cluster", tr_ncluster, tr_id_cluster);
473 sc = m_tuple_track->addIndexedItem("x_cluster", tr_ncluster, tr_x_cluster);
474 sc = m_tuple_track->addIndexedItem("y_cluster", tr_ncluster, tr_y_cluster);
475 sc = m_tuple_track->addIndexedItem("z_cluster", tr_ncluster, tr_z_cluster);
476 sc = m_tuple_track->addIndexedItem("Q_cluster", tr_ncluster, tr_Q_cluster);
477 sc = m_tuple_track->addIndexedItem("phi_cluster", tr_ncluster, tr_phi_cluster);
478 sc = m_tuple_track->addIndexedItem("layer_cluster", tr_ncluster, tr_layer_cluster);
479 sc = m_tuple_track->addIndexedItem("sheet_cluster", tr_ncluster, tr_sheet_cluster);
480 sc = m_tuple_track->addIndexedItem("Is_selected_cluster", tr_ncluster, tr_Is_selected_cluster);
481
482 } else {
483 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple_track) << endmsg;
484 return StatusCode::FAILURE;
485 }
486 }
487
488 IRawDataProviderSvc *irawDataProviderSvc;
489 sc = service("RawDataProviderSvc", irawDataProviderSvc);
490 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc *>(irawDataProviderSvc);
491 if (sc.isFailure()) {
492 log << MSG::FATAL << name() << " Could not load RawDataProviderSvc!" << endreq;
493 return StatusCode::FAILURE;
494 }
495
496 ICgemGeomSvc *icgemGeomSvc;
497 sc = service("CgemGeomSvc", icgemGeomSvc);
498 m_cgemGeomSvc = dynamic_cast<CgemGeomSvc *>(icgemGeomSvc);
499
500 if (sc.isFailure()) {
501 log << MSG::FATAL << "Could not load CgemGeomSvc!" << endreq;
502 return StatusCode::FAILURE;
503 }
504 myNCgemLayers = m_cgemGeomSvc->getNumberOfCgemLayer();
505 for (int i = 0; i < myNCgemLayers; i++) {
506 myCgemLayer[i] = m_cgemGeomSvc->getCgemLayer(i);
507 myNSheets[i] = myCgemLayer[i]->getNumberOfSheet();
508 myRXstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu2();
509 myRVstrips[i] = myCgemLayer[i]->getInnerROfAnodeCu1();
510 }
511 R_layer[0] = (m_cgemGeomSvc->getCgemLayer(0)->getMiddleROfGapD());
512 R_layer[1] = (m_cgemGeomSvc->getCgemLayer(1)->getMiddleROfGapD());
513 R_layer[2] = (m_cgemGeomSvc->getCgemLayer(2)->getMiddleROfGapD());
514 AngleOfStereo[0] = (m_cgemGeomSvc->getCgemLayer(0)->getAngleOfStereo());
515 AngleOfStereo[1] = (m_cgemGeomSvc->getCgemLayer(1)->getAngleOfStereo());
516 AngleOfStereo[2] = (m_cgemGeomSvc->getCgemLayer(2)->getAngleOfStereo());
517 length = m_cgemGeomSvc->getLengthOfCgem();
518 Mp = m_cgemGeomSvc->getMidDriftPtr();
519 Al = m_cgemGeomSvc->getAlignPtr();
520
521 GeoLayer0 = m_cgemGeomSvc->getCgemLayer(0);
522 GeoLayer1 = m_cgemGeomSvc->getCgemLayer(1);
523 GeoLayer2 = m_cgemGeomSvc->getCgemLayer(2);
524
525 pl_00 = m_cgemGeomSvc->getReadoutPlane(0, 0);
526 pl_10 = m_cgemGeomSvc->getReadoutPlane(1, 0);
527 pl_11 = m_cgemGeomSvc->getReadoutPlane(1, 1);
528 pl_20 = m_cgemGeomSvc->getReadoutPlane(2, 0);
529 pl_21 = m_cgemGeomSvc->getReadoutPlane(2, 1);
530
531 cout << "CgemLineFit alignment: is it on? " << Align_Flag << endl;
532 if (Align_Flag) {
533 for (int ilay = 0; ilay < 6; ilay++) {
534 cout << "LAYER " << ilay + 1 << endl;
535 cout << "shift dx " << Al->getDx(ilay) << " dy " << Al->getDy(ilay) << " dz " << Al->getDz(ilay) << endl;
536 cout << "rotation Rx " << Al->getRx(ilay) << " Ry " << Al->getRy(ilay) << " Rz " << Al->getRz(ilay) << endl;
537 cout << "midplane radius " << Mp->getR(int(ilay / 2)) << endl;
538 }
539 }
540
541 ICgemCalibFunSvc *icgemCalibSvc;
542 sc = service("CgemCalibFunSvc", icgemCalibSvc);
543 m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc *>(icgemCalibSvc);
544 if (sc.isFailure()) {
545 log << MSG::FATAL << "Could not load CgemCalibFunSvc!" << endreq;
546 return StatusCode::FAILURE;
547 }
548 sc = service("CgemCalibFunSvc", myCgemCalibSvc);
549 if (sc.isFailure()) {
550 cout << name() << "Could not load MdcCalibFunSvc!" << endl;
551 return sc;
552 }
553 // cout<<"sigma from CgemCalibFunSvc: "<<myCgemCalibSvc->getSigma(0,0,0,0,0,0)<<endl;
554 lutreader = new CgemLUTReader(m_LUTfile);
555 lutreader->ReadLUT();
556
557 return StatusCode::SUCCESS;
558}
559
560// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
561StatusCode CgemLineFit::execute() {
562 MsgStream log(msgSvc(), name());
563 log << MSG::INFO << "in execute()" << endreq;
564 if (debug) cout << "CgemLineFit::execute START" << endl;
565
566 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(), "/Event/EventHeader");
567 if (!eventHeader) {
568 log << MSG::FATAL << "Could not find Event Header" << endreq;
569 return StatusCode::FAILURE;
570 }
571
572 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
573 if (!recCgemClusterCol) {
574 log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
575 return StatusCode::FAILURE;
576 }
577 SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(), "/Event/Digi/CgemDigiCol");
578 if (!cgemDigiCol) {
579 log << MSG::WARNING << "Could not retrieve Cgem digi list" << endreq;
580 return StatusCode::FAILURE;
581 }
582
583 if (debug) cout << "CgemLineFit::execute - read the collection" << endl;
584
585 m_event = eventHeader->eventNumber();
586 m_run = eventHeader->runNumber();
587 event = m_event;
588 run = m_run;
589 tr_run = run;
590 tr_event = event;
591 if ((m_event % 500) == 0) cout << " begin evt : " << event << endl;
592 //cout<<__FILE__<<": evt = "<<event<<endl;
593 //////////////////Collecting Clusters///////////
594
595 // release the memory
596 Vec_layer.clear();
597 Vec_phi.clear();
598 Vec_Z.clear();
599 Vec_v.clear();
600 Vec_layerid.clear();
601 Vec_Virlayerid.clear();
602 Vec_flag.clear();
603 Vec_Q.clear();
604 Vec_sheetid.clear();
605 Vec_m_layer.clear();
606 Vec_m_phi.clear();
607 Vec_m_Z.clear();
608 Vec_m_v.clear();
609 Vec_m_layerid.clear();
610 Vec_m_flag.clear();
611 Vec_m_Q.clear();
612 Vec_m_sheetid.clear();
613 vecclusters.clear();
614 SmartRefVector<RecCgemCluster>().swap(vecclusters);
615
616 _clst_0 = 0;
617 _clst_1 = 0;
618 _clst_2 = 0;
619
620 NXComb = -1;
621 NVComb = -1;
622 // int is_cluster_exsited[6] = {0,0,0,0,0,0};
623 // int nEvent_with_6cluster = 0;
624 // int nEvent_with_4cluster = 0;
625 // for(int i=0;i<6;i++){is_cluster_exsited[i]=0;}
626
627 if (debug) cout << "CgemLineFit::execute - cleaned the variables" << endl;
628
629 TStopwatch gtimer;
630 gtimer.Start();
631 if (debug) cout << "start Straight line fit" << endl;
632 if (Method == 0) {
633 if (!ToyMC()) {
634 tr_Is_fitted = 0;
635 tr_chisq = -1;
636 m_tuple_track->write();
637 if (debug) cout << "CgemLineFit::execute - method 0 failed" << endl;
638 return StatusCode::SUCCESS;
639 }
640 } else if (Method == 1) {
641 if (!Data_Max()) {
642 tr_Is_fitted = 0;
643 tr_chisq = -1;
644 m_tuple_track->write();
645 if (debug) cout << "CgemLineFit::execute - method 1 failed" << endl;
646 return StatusCode::SUCCESS;
647 }
648 } else if (Method == 2) {
649 if (!Data_Closest()) {
650 tr_Is_fitted = 0;
651 tr_chisq = -1;
652 m_tuple_track->write();
653 if (debug) cout << "CgemLineFit::execute - method 2 failed" << endl;
654 return StatusCode::SUCCESS;
655 }
656 } else if (Method == 3) {
657 if (!Loop_All()) {
658 tr_Is_fitted = 0;
659 tr_chisq = -1;
660 m_tuple_track->write();
661 if (debug) cout << "CgemLineFit::execute - method 3 failed" << endl;
662 return StatusCode::SUCCESS;
663 }
664 } else if (Method == 4) {
665 if (!Loop_MaxQ()) {
666 tr_Is_fitted = 0;
667 tr_chisq = -1;
668 m_tuple_track->write();
669 // cout << "Nevent with 6 clusters: " << nEvent_with_6cluster << endl;
670 // cout << "Nevent with 4 clusters: " << nEvent_with_4cluster << endl;
671 if (debug) cout << "CgemLineFit::execute - method 4 failed" << endl;
672 return StatusCode::SUCCESS;
673 }
674 // cout << "Nevent with 6 clusters: " << nEvent_with_6cluster << endl;
675 // cout << "Nevent with 4 clusters: " << nEvent_with_4cluster << endl;
676 }
677
678 if (Vec_layerid.size() < 4 && Run10_Flag) return StatusCode::SUCCESS;
679
680 /////Set initial parameter according to the clusters on the outtermost layer and perform the final fit:
681 /////1, fit phi0, drho; 2, fit tanl z0; 3, fit 4 par. at the same time/////
682
683 if (debug) cout << "CgemLineFit::execute - set initial parameter" << endl;
684
685 TMinuit *fit = Fit(l_outer, fa, 0);
686 double trkpar[4] = {-999, -999, -999, -999};
687 double trkparErr[4] = {-999, -999, -999, -999};
688 Int_t ierflg, npari, nparx, istat3;
689 Double_t fmin, edm, errdef;
690 fit->mnstat(fmin, edm, errdef, npari, nparx, istat3);
691 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
692 double emat[16];
693 fit->mnemat(emat, 4);
694 chisq_3 = fit->fAmin;
695 registerTrack(m_trackList_tds, m_hitList_tds);
696 Store_Trk(fit, 0);
697 delete fit;
698 double chisq_last = chisq_3;
699
700 gtimer.Stop();
701 vector<double> _trkpar = Trans(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
702
703 par0 = _trkpar[0];
704 par1 = _trkpar[1];
705 par2 = _trkpar[2];
706 par3 = _trkpar[3];
707 Err_par0 = trkparErr[0];
708 Err_par1 = trkparErr[1];
709 Err_par2 = trkparErr[2];
710 Err_par3 = trkparErr[3];
711
712 status3 = istat3;
713
714 clst_0 = _clst_0;
715 clst_1 = _clst_1;
716 clst_2 = _clst_2;
717
718 R0 = R_layer[0];
719 R1 = R_layer[1];
720 R2 = R_layer[2];
721
722 // the initial track parameters based on the clusters on the outtermost layer
723 ini_0 = l_outer->dr();
724 ini_1 = l_outer->phi0();
725 ini_2 = l_outer->dz();
726 ini_3 = l_outer->tanl();
727
728 if (debug) cout << "Get_otherIniPar" << endl;
729 // Get the initial track parameters based on the clusters on middle and innermost layers
730 if (TEST_N == 0) Get_OtherIniPar(_clst_0, _clst_1, _clst_2);
731
732 // Store the x,v,z information of clusters on the avialable layers,_clst_0: N clusters on middle layer, _clst_1: N clusters on
733 // innermost layer _clst_0=2; _clst_1=2; _clst_2=0;
734
735 if (debug) cout << "Rec_Clusters" << endl;
736 Rec_Clusters();
737 Rec_ClusterQ();
738 Rec_ClusterSize();
739 fireStripQ();
740 fireStripT();
741 fireStripID();
742 // Store the x,v,z information of clusters by mTPC method on the avialable layers
743 if (Switch_CCmTPC == 1) Rec_Clusters_mTPC();
744
745 if (debug) cout << "Fit_Clusters" << endl;
746 // Store the intersections between fitted line and Cgem
747 Fit_Clusters(trkpar);
748
749 if (debug) cout << "GetIntersection" << endl;
750 // Store the intersections according to the extroplation of the clusters on other layers
751 StraightLine interline(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
752 GetIntersection(&interline);
753
754 if (debug) cout << "GetRes" << endl;
755 // Get the resolution
756 StraightLine l_fit(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
757 GetRes(&l_fit, TEST_N);
758 m_tuple->write();
759
760 tr_drho = trkpar[0];
761 tr_phi0 = trkpar[1];
762 tr_dz0 = trkpar[2];
763 tr_tanLam = trkpar[3];
764 tr_Is_fitted = 1;
765 tr_chisq = chisq_last;
766 m_tuple_track->write();
767
768 delete l_outer;
769 if (debug) cout << "CgemLineFit::execute END" << endl;
770
771 ////////
772 //cout << "NC = " << NC << endl;
773 //SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
774 //if (!recMdcTrackCol) {
775 // log << MSG::WARNING << "Could not retrieve Rec Track list" << endreq;
776 // return StatusCode::FAILURE;
777 //}
778
779 //RecMdcTrackCol::iterator itera = recMdcTrackCol->begin();
780 //cout << "Ncluster = " << (*itera)->getNcluster() << endl;
781 /////////
782 return StatusCode::SUCCESS;
783}
784
785StatusCode CgemLineFit::finalize() {
786 MsgStream log(msgSvc(), name());
787 log << MSG::INFO << "in finalize()" << endreq;
788 return StatusCode::SUCCESS;
789}
790
791//========================================================================================================
792// Method0: for ToyMC study
793//========================================================================================================
794
795bool CgemLineFit::ToyMC() {
796
797 double cl_00(0), cl_01(0), cl_11(0), cl_10(0);
798 int _loop(0);
799 if (!Get_MCInf()) return false;
800 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
801 RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
802 for (; iter != recCgemClusterCol->end(); iter++) {
803 _loop++;
804 RecCgemCluster *cluster = (*iter);
805 int flag = cluster->getflag();
806 if (flag != 2) continue;
807
808 int layerid = cluster->getlayerid();
809 int sheetid = cluster->getsheetid();
810 double _phi = cluster->getrecphi();
811 if (MC) _phi = RealPhi(_phi);
812 double _v = cluster->getrecv();
813 double _Z = cluster->getRecZ();
814 double _Q = cluster->getenergydeposit();
815 int vir_layerid = GetVirLay(layerid, _phi);
816 if (layerid == 2 || ((!Run10_Flag) && layerid == 0 && _phi > acos(-1))) continue;
817
818 Vec_flag.push_back(flag);
819 Vec_layerid.push_back(layerid);
820 Vec_Virlayerid.push_back(vir_layerid);
821 Vec_sheetid.push_back(sheetid);
822 Vec_layer.push_back(R_layer[layerid]);
823
824 Vec_phi.push_back(_phi);
825 Vec_v.push_back(_v);
826 Vec_Z.push_back(_Z);
827
828 Vec_Q.push_back(_Q);
829 cluster->setTrkId(0);
830 vecclusters.push_back(cluster);
831 }
832 int NC = vecclusters.size();
833 if (debug) cout << "NC is " << NC << endl;
834
835 if (Vec_layer.size() < 3) { return false; }
836
837 for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
838 if (Vec_layerid[i_cl] == 2) _clst_2++;
839 if (Vec_layerid[i_cl] == 1) _clst_1++;
840 if (Vec_layerid[i_cl] == 0) _clst_0++;
841 }
842
843 if (debug) cout << "clst1, clst0 are " << _clst_1 << ", " << _clst_0 << endl;
844 OrderClusters();
845 if (_clst_1 == 2) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2); // should be this
846 else return false;
847
848 if (_clst_0 > 2) Filter(0, l_outer);
849
850 return true;
851}
852
853//========================================================================================================
854// Method1: select the cluster with the largest charge on each hemisphere
855//========================================================================================================
856
857bool CgemLineFit::Data_Max() {
858 if (debug) cout << "CgemLineFit::DataMax - START" << endl;
859
860 double cl_00(0), cl_01(0), cl_11(0), cl_10(0), cl_20(0), cl_21(0);
861 int _loop(0);
862 double maxQ_11(0), maxQ_10(0), maxQ_00(0), maxQ_01(0), maxQ_21(0), maxQ_20(0);
863 int max_11(-1), max_10(-1), max_00(-1), max_01(-1), max_21(-1), max_20(-1);
864 double maxv_00(0), maxv_01(0);
865
866 double phi_11(-1), phi_10(-1), phi_00(-1), phi_01(-1), phi_21(-1), phi_20(-1);
867 double z_11(0), z_10(0), z_00(0), z_01(0), z_21(0), z_20(0);
868 double cid_11(-1), cid_10(-1), cid_00(-1), cid_01(-1), cid_21(-1), cid_20(-1);
869 double Xid_11(-1), Xid_10(-1), Xid_00(-1), Xid_01(-1), Xid_21(-1), Xid_20(-1);
870 double Vid_11(-1), Vid_10(-1), Vid_00(-1), Vid_01(-1), Vid_21(-1), Vid_20(-1);
871
872 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
873 RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
874
875 int ic(0);
876 int ic_tot(0);
877 if (debug) cout << "CgemLineFit::DataMax - First loop on cluster" << endl;
878 for (; iter != recCgemClusterCol->end(); iter++) {
879 _loop++;
880 if (debug) cout << _loop << endl;
881 RecCgemCluster *cluster = (*iter);
882 int flag = cluster->getflag();
883 int layerid = cluster->getlayerid();
884 double _Q = cluster->getenergydeposit();
885 double _phi = cluster->getrecphi();
886 double _v = cluster->getrecv();
887 int sheetid = cluster->getsheetid();
888
889 int clusterid = cluster->getclusterid();
890 int clusterflagb = cluster->getclusterflagb();
891 int clusterflage = cluster->getclusterflage();
892 double x = R_layer[layerid] * cos(_phi);
893 double y = R_layer[layerid] * sin(_phi);
894 double z = cluster->getRecZ();
895 if (debug) cout << " flag: " << flag << endl;
896 if (debug)
897 cout << "CgemLineFit::DataMax - cluster info - ID: " << clusterid << " layer: " << layerid << " sheet: " << sheetid
898 << " Xid: " << clusterflagb << " Vid: " << clusterflage << " Phi: " << _phi << " Z: " << z << " X: " << x
899 << " Y: " << y << " Q: " << _Q << endl;
900
901 if (_loop > MAX_Clusters) break;
902 if (flag != 2) continue; // only cluster 2D
903 if (_Q < MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
904
905 tr_id_cluster[ic] = clusterid;
906 tr_x_cluster[ic] = x;
907 tr_y_cluster[ic] = y;
908 tr_z_cluster[ic] = z;
909 tr_Q_cluster[ic] = _Q;
910 // tr_phi_cluster[ic] = _phi;
911 // tr_layer_cluster[ic] = layerid;
912 // tr_sheet_cluster[ic] = sheetid;
913
914 ic_tot++;
915 ic++;
916
917 if (layerid == 2 && sheetid == 1 && maxQ_21 < _Q) {
918 maxQ_21 = _Q;
919 max_21 = _loop;
920 phi_21 = _phi;
921 z_21 = z;
922 cid_21 = clusterid;
923 Xid_21 = clusterflagb;
924 Vid_21 = clusterflage;
925 }
926 if (layerid == 2 && sheetid == 0 && maxQ_20 < _Q) {
927 maxQ_20 = _Q;
928 max_20 = _loop;
929 phi_20 = _phi;
930 z_20 = z;
931 cid_20 = clusterid;
932 Xid_20 = clusterflagb;
933 Vid_20 = clusterflage;
934 }
935 if (layerid == 1 && sheetid == 1 && maxQ_11 < _Q) {
936 maxQ_11 = _Q;
937 max_11 = _loop;
938 phi_11 = _phi;
939 z_11 = z;
940 cid_11 = clusterid;
941 Xid_11 = clusterflagb;
942 Vid_11 = clusterflage;
943 }
944 if (layerid == 1 && sheetid == 0 && maxQ_10 < _Q) {
945 maxQ_10 = _Q;
946 max_10 = _loop;
947 phi_10 = _phi;
948 z_10 = z;
949 cid_10 = clusterid;
950 Xid_10 = clusterflagb;
951 Vid_10 = clusterflage;
952 }
953 if (layerid == 0 && sheetid == 0 && _phi > 0 && maxQ_01 < _Q && maxv_00 != _v) {
954 maxQ_01 = _Q;
955 max_01 = _loop;
956 phi_01 = _phi;
957 z_01 = z;
958 cid_01 = clusterid;
959 Xid_01 = clusterflagb;
960 Vid_01 = clusterflage;
961 }
962 if (layerid == 0 && sheetid == 0 && _phi < 0 && maxQ_00 < _Q && maxv_00 != _v) {
963 maxQ_00 = _Q;
964 max_00 = _loop;
965 phi_00 = _phi;
966 z_00 = z;
967 cid_00 = clusterid;
968 Xid_00 = clusterflagb;
969 Vid_00 = clusterflage;
970 }
971
972 if (layerid == 2 && sheetid == 1) cl_21++;
973 if (layerid == 2 && sheetid == 0) cl_20++;
974 if (layerid == 1 && sheetid == 1) cl_11++;
975 if (layerid == 1 && sheetid == 0) cl_10++;
976 if (layerid == 0 && sheetid == 0 && _phi > 0) cl_01++;
977 if (layerid == 0 && sheetid == 0 && _phi < 0) cl_00++;
978 }
979
980 if (debug) cout << "CgemLineFit::DataMax - First loop on cluster ended - number of cluster: " << ic << endl;
981
982 tr_ncluster = ic;
983 for (int i = 0; i < ic; i++) {
984 tr_Is_selected_cluster[i] = 0;
985 if (i == (max_00 - 1) && maxQ_00 > 0) tr_Is_selected_cluster[i] = 1;
986 if (i == (max_01 - 1) && maxQ_01 > 0) tr_Is_selected_cluster[i] = 1;
987 if (i == (max_10 - 1) && maxQ_10 > 0) tr_Is_selected_cluster[i] = 1;
988 if (i == (max_11 - 1) && maxQ_11 > 0) tr_Is_selected_cluster[i] = 1;
989 if (i == (max_20 - 1) && maxQ_20 > 0) tr_Is_selected_cluster[i] = 1;
990 if (i == (max_21 - 1) && maxQ_21 > 0) tr_Is_selected_cluster[i] = 1;
991 }
992
993 clst_21 = cl_21;
994 clst_20 = cl_20;
995 clst_11 = cl_11;
996 clst_01 = cl_01;
997 clst_10 = cl_10;
998 clst_00 = cl_00;
999 _loop = 0;
1000
1001 if (debug)
1002 cout << "CgemLineFit::DataMax - number of cluster per sheet/layer: " << clst_00 << " " << clst_10 << " " << clst_01 << " "
1003 << clst_11 << " " << clst_20 << " " << clst_21 << endl;
1004
1005 iter = recCgemClusterCol->begin();
1006 for (; iter != recCgemClusterCol->end(); iter++) {
1007 _loop++;
1008 RecCgemCluster *cluster = (*iter);
1009 int flag = cluster->getflag();
1010 if (flag != 2) continue;
1011
1012 int layerid = cluster->getlayerid();
1013 int sheetid = cluster->getsheetid();
1014 double _phi = cluster->getrecphi();
1015 double _v = cluster->getrecv();
1016 double _Z = cluster->getRecZ();
1017 double _Q = cluster->getenergydeposit();
1018 double t_phi = cluster->getrecphi_mTPC();
1019 double t_v = cluster->getrecv_mTPC();
1020 double t_Z = cluster->getRecZ_mTPC();
1021 int vir_layerid = GetVirLay(layerid, _phi);
1022
1023 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
1024 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
1025 double merged_phi = cluster_x->get_merge_phi();
1026 RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
1027 double merged_v = cluster_v->get_merge_v();
1028
1029 // Added to be checked - farinelli 20240106
1030 int clusterflagb = cluster->getclusterflagb();
1031 int clusterflage = cluster->getclusterflage();
1032 RecCgemClusterCol::iterator id_XCluster = recCgemClusterCol->begin() + clusterflagb;
1033 RecCgemCluster *xcluster = (*id_XCluster);
1034 int xclustersize = abs(xcluster->getclusterflage() - xcluster->getclusterflagb()) + 1;
1035 double X_Q = xcluster->getenergydeposit();
1036 RecCgemClusterCol::iterator id_VCluster = recCgemClusterCol->begin() + clusterflage;
1037 RecCgemCluster *vcluster = (*id_VCluster);
1038 int vclustersize = abs(vcluster->getclusterflage() - vcluster->getclusterflagb()) + 1;
1039 double V_Q = vcluster->getenergydeposit();
1040 int xclusterb = xcluster->getclusterflagb();
1041 int xclustere = xcluster->getclusterflage();
1042 int vclusterb = vcluster->getclusterflagb();
1043 int vclustere = vcluster->getclusterflage();
1044 int layerid_str, sheetid_str, stripid_str, time_chnnel;
1045 int nxstrips = 0;
1046 int nvstrips = 0;
1047 double energydeposit;
1048 double Q_fC, T_ns, Tf_ns;
1049 SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(), "/Event/Digi/CgemDigiCol");
1050 CgemDigiCol::iterator iter_digi = cgemDigiCol->begin();
1051 Vec_x_stripsQ.clear();
1052 Vec_v_stripsQ.clear();
1053 Vec_x_stripsT.clear();
1054 Vec_v_stripsT.clear();
1055 Vec_x_stripsTf.clear();
1056 Vec_v_stripsTf.clear();
1057 Vec_x_stripsid.clear();
1058 Vec_v_stripsid.clear();
1059 for (; iter_digi != cgemDigiCol->end(); iter_digi++) {
1060 layerid_str = CgemID::layer((*iter_digi)->identify());
1061 sheetid_str = CgemID::sheet((*iter_digi)->identify());
1062 stripid_str = CgemID::strip((*iter_digi)->identify());
1063 bool striptype = CgemID::is_xstrip((*iter_digi)->identify());
1064 int flagxv;
1065 if (striptype == true) flagxv = 0;
1066 else flagxv = 1;
1067 if (layerid == layerid_str && sheetid == sheetid_str) {
1068 Q_fC = (*iter_digi)->getCharge_fc();
1069 T_ns = (*iter_digi)->getTime_ns();
1070 Tf_ns = get_Time(iter_digi);
1071 if (flagxv == 0) {
1072 if (stripid_str >= xclusterb && stripid_str <= xclustere) {
1073 if (xclustersize < 100) {
1074 Vec_x_stripsid.push_back(stripid_str);
1075 Vec_x_stripsQ.push_back(Q_fC);
1076 Vec_x_stripsT.push_back(T_ns);
1077 Vec_x_stripsTf.push_back(Tf_ns);
1078 nxstrips++;
1079 }
1080 }
1081 }
1082 if (flagxv == 1) {
1083 if (stripid_str >= vclusterb && stripid_str <= vclustere) {
1084 if (vclustersize < 100) {
1085 Vec_v_stripsid.push_back(stripid_str);
1086 Vec_v_stripsQ.push_back(Q_fC);
1087 Vec_v_stripsT.push_back(T_ns);
1088 Vec_v_stripsTf.push_back(Tf_ns);
1089 nvstrips++;
1090 }
1091 }
1092 }
1093 }
1094 }
1095 // End to be checked
1096
1097 if (debug) cout << "cl ID: " << _loop << endl;
1098 if (!(_loop == max_00 || _loop == max_10 || _loop == max_01 || _loop == max_11 || _loop == max_20 || _loop == max_21))
1099 continue;
1100 if (debug) cout << "is selected" << endl;
1101
1102 Vec_flag.push_back(flag);
1103 Vec_layerid.push_back(layerid);
1104 Vec_sheetid.push_back(sheetid);
1105 Vec_layer.push_back(R_layer[layerid]);
1106 Vec_Virlayerid.push_back(vir_layerid);
1107
1108 Vec_m_phi.push_back(t_phi);
1109 Vec_m_v.push_back(t_v);
1110 Vec_m_Z.push_back(t_Z);
1111
1112 Vec_xCluSize.push_back(xclustersize);
1113 Vec_vCluSize.push_back(vclustersize);
1114 Vec_XQ.push_back(X_Q);
1115 Vec_VQ.push_back(V_Q);
1116 Vec_XstripQ.push_back(Vec_x_stripsQ);
1117 Vec_VstripQ.push_back(Vec_v_stripsQ);
1118 Vec_XstripT.push_back(Vec_x_stripsT);
1119 Vec_VstripT.push_back(Vec_v_stripsT);
1120 Vec_XstripTf.push_back(Vec_x_stripsTf);
1121 Vec_VstripTf.push_back(Vec_v_stripsTf);
1122 Vec_XstripID.push_back(Vec_x_stripsid);
1123 Vec_VstripID.push_back(Vec_v_stripsid);
1124
1125 if (Switch_CCmTPC == 0) {
1126 Vec_phi.push_back(_phi);
1127 Vec_v.push_back(_v);
1128 Vec_Z.push_back(_Z);
1129
1130 } else if (Switch_CCmTPC == 1) {
1131 Vec_phi.push_back(t_phi);
1132 Vec_v.push_back(t_v);
1133 Vec_Z.push_back(t_Z);
1134
1135 } else if (Switch_CCmTPC == 2) {
1136
1137 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
1138 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
1139 Vec_phi.push_back(merged_phi);
1140 Vec_v.push_back(merged_v);
1141 Vec_Z.push_back(merged_Z);
1142 }
1143
1144 Vec_Q.push_back(_Q);
1145 cluster->setTrkId(0);
1146 vecclusters.push_back(cluster);
1147 }
1148 NC = vecclusters.size();
1149
1150 if (debug) cout << "CgemLineFit::DataMax - Number of points: " << Vec_layer.size() << endl;
1151
1152 // if(Vec_layer.size()!=6) return false;
1153 if (Vec_layer.size() < 2) return false;
1154 NXComb = 1;
1155 NVComb = 1;
1156
1157 for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
1158 if (Vec_layerid[i_cl] == 2) _clst_2++;
1159 if (Vec_layerid[i_cl] == 1) _clst_1++;
1160 if (Vec_layerid[i_cl] == 0) _clst_0++;
1161 }
1162
1163 OrderClusters();
1164
1165 // if(debug) cout<<"CgemLineFit::DataMax - _clst_2: "<<_clst_2<<endl;
1166 // if(_clst_1==2) l_outer=IniPar(Vec_phi[0],Vec_Z[0],3,Vec_phi[1],Vec_Z[1],2);
1167 // changed to be checked - farinelli 20240106
1168 if (debug) cout << "CgemLineFit::DataMax - Vec_layerid.size(): " << Vec_layerid.size() << endl;
1169 if (Vec_layerid.size() >= 2)
1170 l_outer = IniPar(Vec_phi[0], Vec_Z[0], 2 * Vec_layerid[0], Vec_phi[1], Vec_Z[1], 2 * Vec_layerid[1]);
1171 else return false;
1172
1173 return true;
1174}
1175
1176//========================================================================================================
1177// Method2: select the cluster closest to the intersection. probably not be finished?
1178//========================================================================================================
1179
1180bool CgemLineFit::Data_Closest() {
1181 double cl_00(0), cl_01(0), cl_11(0), cl_10(0);
1182 int _loop(0);
1183 double maxQ_11(0), maxQ_10(0), maxQ_00(0);
1184 int max_11(-1), max_10(-1), max_00(-1);
1185 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
1186 RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
1187
1188 for (; iter != recCgemClusterCol->end(); iter++) {
1189 _loop++;
1190 RecCgemCluster *cluster = (*iter);
1191 int flag = cluster->getflag();
1192 int layerid = cluster->getlayerid();
1193 double _Q = cluster->getenergydeposit();
1194 int sheetid = cluster->getsheetid();
1195 double _phi = cluster->getrecphi();
1196 if (flag != 2) continue;
1197 if (layerid == 1 && sheetid == 1 && maxQ_11 < _Q) {
1198 maxQ_11 = _Q;
1199 max_11 = _loop;
1200 }
1201 if (layerid == 1 && sheetid == 0 && maxQ_10 < _Q) {
1202 maxQ_10 = _Q;
1203 max_10 = _loop;
1204 }
1205 if (layerid == 1 && sheetid == 1) cl_11++;
1206 if (layerid == 0 && sheetid == 0 && _phi > 0) cl_01++;
1207 if (layerid == 0 && sheetid == 0 && _phi < 0) cl_00++;
1208 if (layerid == 1 && sheetid == 0) cl_10++;
1209 }
1210 clst_11 = cl_11;
1211 clst_01 = cl_01;
1212 clst_10 = cl_10;
1213 clst_00 = cl_00;
1214 _loop = 0;
1215 iter = recCgemClusterCol->begin();
1216 for (; iter != recCgemClusterCol->end(); iter++) {
1217 _loop++;
1218 RecCgemCluster *cluster = (*iter);
1219 int flag = cluster->getflag();
1220 if (flag != 2) continue;
1221
1222 int layerid = cluster->getlayerid();
1223 int sheetid = cluster->getsheetid();
1224 double _phi = cluster->getrecphi();
1225 double _v = cluster->getrecv();
1226 double _Z = cluster->getRecZ();
1227 double _Q = cluster->getenergydeposit();
1228 double t_phi = cluster->getrecphi_mTPC();
1229 double t_v = cluster->getrecv_mTPC();
1230 double t_Z = cluster->getRecZ_mTPC();
1231 int vir_layerid = GetVirLay(layerid, _phi); // Added by Farinelli Riccardo in 2023.12.29
1232
1233 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
1234 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
1235 double merged_phi = cluster_x->get_merge_phi();
1236 RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
1237 double merged_v = cluster_v->get_merge_v();
1238
1239 if (_loop != max_11 && _loop != max_10 && layerid != 0) continue;
1240 Vec_flag.push_back(flag);
1241 Vec_layerid.push_back(layerid);
1242 Vec_sheetid.push_back(sheetid);
1243 Vec_layer.push_back(R_layer[layerid]);
1244 Vec_Virlayerid.push_back(vir_layerid); // Added by Farinelli Riccardo in 2023.12.29
1245
1246 Vec_m_phi.push_back(t_phi);
1247 Vec_m_v.push_back(t_v);
1248 Vec_m_Z.push_back(t_Z);
1249
1250 if (Switch_CCmTPC == 0) {
1251 Vec_phi.push_back(_phi);
1252 Vec_v.push_back(_v);
1253 Vec_Z.push_back(_Z);
1254 } else if (Switch_CCmTPC == 1) {
1255 Vec_phi.push_back(t_phi);
1256 Vec_v.push_back(t_v);
1257 Vec_Z.push_back(t_Z);
1258 } else if (Switch_CCmTPC == 2) {
1259
1260 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
1261 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
1262 Vec_phi.push_back(merged_phi);
1263 Vec_v.push_back(merged_v);
1264 Vec_Z.push_back(merged_Z);
1265 }
1266
1267 Vec_Q.push_back(_Q);
1268 cluster->setTrkId(0);
1269 vecclusters.push_back(cluster);
1270 }
1271 NC = vecclusters.size();
1272
1273 if (Vec_layer.size() < 3) { return false; }
1274 for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
1275 if (Vec_layerid[i_cl] == 2) _clst_2++;
1276 if (Vec_layerid[i_cl] == 1) _clst_1++;
1277 if (Vec_layerid[i_cl] == 0) _clst_0++;
1278 }
1279
1280 OrderClusters();
1281
1282 if (_clst_1 == 2) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 1, Vec_phi[1], Vec_Z[1], 1);
1283 else return false;
1284
1285 if (!Layer_cross(0, l_outer)) return false;
1286 if (_clst_0 > 1 || (Run10_Flag && _clst_0 > 2)) Filter(0, l_outer);
1287
1288 return true;
1289}
1290
1291//========================================================================================================
1292// Method3: loop all the clusters then select the combination giving the smallest chisq
1293//========================================================================================================
1294
1295bool CgemLineFit::Loop_All() {
1296 double cl_00(0), cl_01(0), cl_11(0), cl_10(0);
1297 int _loop(0);
1298 double maxQ_11(0), maxQ_10(0), maxQ_00(0);
1299 int max_11(-1), max_10(-1), max_00(-1);
1300 double C_00(0), C_10(0), C_11(0), C_01(0);
1301 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; //
1302 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; //
1303 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; //
1304 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; //
1305
1306 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
1307 RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
1308 vector<int> L_00, L_01, L_10, L_11;
1309 int ic(0);
1310 for (; iter != recCgemClusterCol->end(); iter++) {
1311 _loop++;
1312 RecCgemCluster *cluster = (*iter);
1313
1314 int flag = cluster->getflag();
1315 int layerid = cluster->getlayerid();
1316 double _Q = cluster->getenergydeposit();
1317 double _phi = cluster->getrecphi();
1318 double _v = cluster->getrecv();
1319 int sheetid = cluster->getsheetid();
1320 int clusterid = cluster->getclusterid();
1321 int clusterflagb = cluster->getclusterflagb();
1322 int clusterflage = cluster->getclusterflage();
1323 double x = R_layer[layerid] * cos(_phi);
1324 double y = R_layer[layerid] * sin(_phi);
1325 double z = cluster->getRecZ();
1326 if (_loop > MAX_Clusters) break;
1327 ic++;
1328 if (flag != 0) continue;
1329
1330 if (layerid == 1 && sheetid == 1 && TEST_N != 1) { L_11.push_back(_loop); }
1331 if (layerid == 1 && sheetid == 0 && TEST_N != 4) { L_10.push_back(_loop); }
1332 if (layerid == 0 && sheetid == 0 && _phi > 0 && TEST_N != 2) { L_01.push_back(_loop); }
1333 if (layerid == 0 && sheetid == 0 && _phi < 0 && TEST_N != 3) { L_00.push_back(_loop); }
1334 }
1335
1336 if (TEST_N == 1) L_11.push_back(-1);
1337 if (TEST_N == 2) L_01.push_back(-1);
1338 if (TEST_N == 3) L_00.push_back(-1);
1339 if (TEST_N == 4) L_10.push_back(-1);
1340
1341 double Min_chi(1E10);
1342
1343 if ((L_00.size() * L_01.size() * L_10.size() * L_11.size()) > MAX_COMB) { return false; }
1344
1345 NXComb = L_00.size() * L_01.size() * L_10.size() * L_11.size();
1346
1347 for (int i_11 = 0; i_11 < L_11.size(); i_11++) {
1348 for (int i_10 = 0; i_10 < L_10.size(); i_10++) {
1349 for (int i_00 = 0; i_00 < L_00.size(); i_00++) {
1350 for (int i_01 = 0; i_01 < L_01.size(); i_01++) {
1351
1352 Vec_layer.clear();
1353 Vec_phi.clear();
1354 Vec_Z.clear();
1355 Vec_v.clear();
1356 Vec_layerid.clear();
1357 Vec_Virlayerid.clear();
1358 Vec_flag.clear();
1359 Vec_Q.clear();
1360 Vec_sheetid.clear();
1361 Vec_m_layer.clear();
1362 Vec_m_phi.clear();
1363 Vec_m_Z.clear();
1364 Vec_m_v.clear();
1365 Vec_m_layerid.clear();
1366 Vec_m_flag.clear();
1367 Vec_m_Q.clear();
1368 Vec_m_sheetid.clear();
1369 _loop = 0;
1370 iter = recCgemClusterCol->begin();
1371 for (; iter != recCgemClusterCol->end(); iter++) {
1372 _loop++;
1373 RecCgemCluster *cluster = (*iter);
1374 int flag = cluster->getflag();
1375 if (!(_loop == L_11[i_11] || _loop == L_10[i_10] || _loop == L_00[i_00] || _loop == L_01[i_01])) continue;
1376 int layerid = cluster->getlayerid();
1377 int sheetid = cluster->getsheetid();
1378 double _phi = cluster->getrecphi();
1379 double _v = cluster->getrecv();
1380 double _Z = cluster->getRecZ();
1381 double _Q = cluster->getenergydeposit();
1382 double t_phi = cluster->getrecphi_mTPC();
1383 double t_v = cluster->getrecv_mTPC();
1384 double t_Z = cluster->getRecZ_mTPC();
1385
1386 double merged_phi = cluster->get_merge_phi();
1387 double merged_v = cluster->get_merge_v();
1388
1389 int vir_layerid = GetVirLay(layerid, _phi);
1390
1391 Vec_flag.push_back(flag);
1392 Vec_layerid.push_back(layerid);
1393 Vec_Virlayerid.push_back(vir_layerid);
1394 Vec_sheetid.push_back(sheetid);
1395 Vec_layer.push_back(R_layer[layerid]);
1396
1397 Vec_m_phi.push_back(t_phi);
1398 Vec_m_v.push_back(t_v);
1399 Vec_m_Z.push_back(t_Z);
1400
1401 if (Switch_CCmTPC == 0) {
1402 Vec_phi.push_back(_phi);
1403 Vec_v.push_back(_v);
1404 Vec_Z.push_back(_Z);
1405
1406 } else if (Switch_CCmTPC == 1) {
1407 Vec_phi.push_back(t_phi);
1408 Vec_v.push_back(t_v);
1409 Vec_Z.push_back(t_Z);
1410
1411 } else if (Switch_CCmTPC == 2) {
1412
1413 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
1414 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
1415 Vec_phi.push_back(merged_phi);
1416 Vec_v.push_back(merged_v);
1417 Vec_Z.push_back(merged_Z);
1418 }
1419
1420 } // one combo
1421 if (Vec_layerid.size() != 4 && TEST_N == 0) continue;
1422 else if (Vec_layerid.size() != 3 && TEST_N > 0) continue;
1423
1424 OrderClusters();
1425
1426 StraightLine *l_outer_tmp;
1427
1428 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);
1429 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);
1430 else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
1431
1432 TMinuit *_fit = Fit(l_outer_tmp, fa, 1);
1433 double _chi = _fit->fAmin;
1434 if (_chi < Min_chi) {
1435 Min_chi = _chi;
1436 if (TEST_N == 0) C_00 = Vec_phi[3];
1437 else C_00 = -99;
1438 C_01 = Vec_phi[2];
1439 C_10 = Vec_phi[1];
1440 C_11 = Vec_phi[0];
1441 }
1442 delete _fit;
1443 delete l_outer_tmp;
1444 }
1445 }
1446 }
1447 }
1448
1449 iter = recCgemClusterCol->begin();
1450 L_00.clear();
1451 L_01.clear();
1452 L_10.clear();
1453 L_11.clear();
1454 _loop = 0;
1455 int ic_tot(0);
1456 for (; iter != recCgemClusterCol->end(); iter++) {
1457 _loop++;
1458 RecCgemCluster *cluster = (*iter);
1459 int flag = cluster->getflag();
1460 if (flag != 2) continue;
1461 ic_tot++;
1462 int layerid = cluster->getlayerid();
1463 double _Q = cluster->getenergydeposit();
1464
1465 double _phi = 99;
1466 if (Switch_CCmTPC == 0) _phi = cluster->getrecphi();
1467 else if (Switch_CCmTPC == 1) _phi = cluster->getrecphi_mTPC();
1468 else if (Switch_CCmTPC == 2) {
1469 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
1470 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
1471 _phi = cluster_x->get_merge_phi();
1472 }
1473
1474 int sheetid = cluster->getsheetid();
1475 if (_Q < MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
1476 if (C_00 == _phi || C_01 == _phi || C_11 == _phi || C_10 == _phi) {
1477 if (layerid == 1 && sheetid == 1) { L_11.push_back(_loop); }
1478 if (layerid == 1 && sheetid == 0) { L_10.push_back(_loop); }
1479 if (layerid == 0 && sheetid == 0 && _phi > 0) { L_01.push_back(_loop); }
1480 if (layerid == 0 && sheetid == 0 && _phi < 0) { L_00.push_back(_loop); }
1481 }
1482 }
1483 Min_chi = 1E10;
1484
1485 if (TEST_N == 1) L_11.push_back(-1);
1486 if (TEST_N == 2) L_01.push_back(-1);
1487 if (TEST_N == 3) L_00.push_back(-1);
1488 if (TEST_N == 4) L_10.push_back(-1);
1489
1490 clst_11 = L_11.size();
1491 clst_01 = L_01.size();
1492 clst_10 = L_10.size();
1493 clst_00 = L_00.size();
1494 if ((L_00.size() * L_01.size() * L_10.size() * L_11.size()) > MAX_COMB) {
1495 cout << "xv size much " << endl;
1496 return false;
1497 }
1498
1499 NVComb = L_00.size() * L_01.size() * L_10.size() * L_11.size();
1500
1501 for (int i_11 = 0; i_11 < L_11.size(); i_11++) {
1502 for (int i_10 = 0; i_10 < L_10.size(); i_10++) {
1503 for (int i_00 = 0; i_00 < L_00.size(); i_00++) {
1504 for (int i_01 = 0; i_01 < L_01.size(); i_01++) {
1505
1506 Vec_layer.clear();
1507 Vec_phi.clear();
1508 Vec_Z.clear();
1509 Vec_v.clear();
1510 Vec_layerid.clear();
1511 Vec_Virlayerid.clear();
1512 Vec_flag.clear();
1513 Vec_Q.clear();
1514 Vec_sheetid.clear();
1515 Vec_m_layer.clear();
1516 Vec_m_phi.clear();
1517 Vec_m_Z.clear();
1518 Vec_m_v.clear();
1519 Vec_m_layerid.clear();
1520 Vec_m_flag.clear();
1521 Vec_m_Q.clear();
1522 Vec_m_sheetid.clear();
1523 _loop = 0;
1524 iter = recCgemClusterCol->begin();
1525 for (; iter != recCgemClusterCol->end(); iter++) {
1526 _loop++;
1527 RecCgemCluster *cluster = (*iter);
1528 int flag = cluster->getflag();
1529 if (flag != 2) continue;
1530
1531 int layerid = cluster->getlayerid();
1532 int sheetid = cluster->getsheetid();
1533 double _phi = cluster->getrecphi();
1534 double _v = cluster->getrecv();
1535 double _Z = cluster->getRecZ();
1536 double _Q = cluster->getenergydeposit();
1537 double t_phi = cluster->getrecphi_mTPC();
1538 double t_v = cluster->getrecv_mTPC();
1539 double t_Z = cluster->getRecZ_mTPC();
1540
1541 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
1542 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
1543 double merged_phi = cluster_x->get_merge_phi();
1544 RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
1545 double merged_v = cluster_v->get_merge_v();
1546
1547 int vir_layerid = GetVirLay(layerid, _phi);
1548
1549 if (!(_loop == L_11[i_11] || _loop == L_10[i_10] || _loop == L_00[i_00] || _loop == L_01[i_01])) continue;
1550
1551 Vec_flag.push_back(flag);
1552 Vec_layerid.push_back(layerid);
1553 Vec_Virlayerid.push_back(vir_layerid);
1554 Vec_sheetid.push_back(sheetid);
1555 Vec_layer.push_back(R_layer[layerid]);
1556
1557 Vec_m_phi.push_back(t_phi);
1558 Vec_m_v.push_back(t_v);
1559 Vec_m_Z.push_back(t_Z);
1560
1561 if (Switch_CCmTPC == 0) {
1562 Vec_phi.push_back(_phi);
1563 Vec_v.push_back(_v);
1564 Vec_Z.push_back(_Z);
1565
1566 } else if (Switch_CCmTPC == 1) {
1567 Vec_phi.push_back(t_phi);
1568 Vec_v.push_back(t_v);
1569 Vec_Z.push_back(t_Z);
1570
1571 } else if (Switch_CCmTPC == 2) {
1572
1573 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
1574 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
1575 Vec_phi.push_back(merged_phi);
1576 Vec_v.push_back(merged_v);
1577 Vec_Z.push_back(merged_Z);
1578 }
1579
1580 } // one combo
1581
1582 if (Vec_layerid.size() != 4 && TEST_N == 0) continue;
1583 else if (Vec_layerid.size() != 3 && TEST_N > 0) continue;
1584
1585 OrderClusters();
1586
1587 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] ||
1588 Vec_v[1] == Vec_v[3] || Vec_v[2] == Vec_v[3])
1589 continue;
1590 StraightLine *l_outer_tmp;
1591
1592 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);
1593 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);
1594 else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
1595
1596 TMinuit *_fit = Fit(l_outer_tmp, fa, 2);
1597 double _chi = _fit->fAmin;
1598 if (_chi < Min_chi) {
1599 Min_chi = _chi;
1600 C_00 = L_00[i_00];
1601 C_01 = L_01[i_01];
1602 C_10 = L_10[i_10];
1603 C_11 = L_11[i_11];
1604 }
1605 delete _fit;
1606 delete l_outer_tmp;
1607 }
1608 }
1609 }
1610 }
1611
1612 Vec_layer.clear();
1613 Vec_phi.clear();
1614 Vec_Z.clear();
1615 Vec_v.clear();
1616 Vec_layerid.clear();
1617 Vec_Virlayerid.clear();
1618 Vec_flag.clear();
1619 Vec_Q.clear();
1620 Vec_sheetid.clear();
1621 Vec_m_layer.clear();
1622 Vec_m_phi.clear();
1623 Vec_m_Z.clear();
1624 Vec_m_v.clear();
1625 Vec_m_layerid.clear();
1626 Vec_m_flag.clear();
1627 Vec_m_Q.clear();
1628 Vec_m_sheetid.clear();
1629 _loop = 0;
1630 ic = 0;
1631 iter = recCgemClusterCol->begin();
1632 for (; iter != recCgemClusterCol->end(); iter++) {
1633 _loop++;
1634 RecCgemCluster *cluster = (*iter);
1635 int flag = cluster->getflag();
1636 int layerid = cluster->getlayerid();
1637 int sheetid = cluster->getsheetid();
1638 double _phi = cluster->getrecphi();
1639 double _v = cluster->getrecv();
1640 double _Z = cluster->getRecZ();
1641 double _Q = cluster->getenergydeposit();
1642 double t_phi = cluster->getrecphi_mTPC();
1643 double t_v = cluster->getrecv_mTPC();
1644 double t_Z = cluster->getRecZ_mTPC();
1645 int clusterid = cluster->getclusterid();
1646 int clusterflagb = cluster->getclusterflagb();
1647 int clusterflage = cluster->getclusterflage();
1648 double x = R_layer[layerid] * cos(_phi);
1649 double y = R_layer[layerid] * sin(_phi);
1650 double z = cluster->getRecZ();
1651 int vir_layerid = GetVirLay(layerid, _phi);
1652
1653 if (_loop > MAX_Clusters) break;
1654 if (flag != 2) continue;
1655
1656 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
1657 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
1658 double merged_phi = cluster_x->get_merge_phi();
1659 RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
1660 double merged_v = cluster_v->get_merge_v();
1661
1662 tr_id_cluster[ic] = clusterid;
1663 tr_x_cluster[ic] = x;
1664 tr_y_cluster[ic] = y;
1665 tr_z_cluster[ic] = z;
1666 tr_Q_cluster[ic] = _Q;
1667 // tr_phi_cluster[ic] = _phi;
1668 // tr_layer_cluster[ic] = layerid;
1669 // tr_sheet_cluster[ic] = sheetid;
1670 tr_Is_selected_cluster[ic] = 0;
1671 ic++;
1672 if (!(_loop == C_00 || _loop == C_10 || _loop == C_01 || _loop == C_11)) continue;
1673
1674 tr_Is_selected_cluster[ic - 1] = 1;
1675
1676 Vec_flag.push_back(flag);
1677 Vec_layerid.push_back(layerid);
1678 Vec_Virlayerid.push_back(vir_layerid);
1679 Vec_sheetid.push_back(sheetid);
1680 Vec_layer.push_back(R_layer[layerid]);
1681
1682 Vec_m_phi.push_back(t_phi);
1683 Vec_m_v.push_back(t_v);
1684 Vec_m_Z.push_back(t_Z);
1685
1686 if (Switch_CCmTPC == 0) {
1687 Vec_phi.push_back(_phi);
1688 Vec_v.push_back(_v);
1689 Vec_Z.push_back(_Z);
1690
1691 } else if (Switch_CCmTPC == 1) {
1692 Vec_phi.push_back(t_phi);
1693 Vec_v.push_back(t_v);
1694 Vec_Z.push_back(t_Z);
1695
1696 } else if (Switch_CCmTPC == 2) {
1697
1698 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
1699 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
1700 Vec_phi.push_back(merged_phi);
1701 Vec_v.push_back(merged_v);
1702 Vec_Z.push_back(merged_Z);
1703 }
1704
1705 cluster->setTrkId(0);
1706 vecclusters.push_back(cluster);
1707
1708 } // one combo
1709 NC = vecclusters.size();
1710 tr_ncluster = ic;
1711
1712 for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
1713 if (Vec_layerid[i_cl] == 2) _clst_2++;
1714 if (Vec_layerid[i_cl] == 1) _clst_1++;
1715 if (Vec_layerid[i_cl] == 0) _clst_0++;
1716 }
1717
1718 if (Vec_layerid.size() != 4 && TEST_N == 0) return false;
1719 else if (Vec_layerid.size() != 3 && TEST_N > 0) return false;
1720
1721 OrderClusters();
1722
1723 if (TEST_N == 1 || TEST_N == 4) l_outer = IniPar(Vec_phi[1], Vec_Z[1], 1, Vec_phi[2], Vec_Z[2], 0);
1724 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);
1725 else if (TEST_N == 0) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1], 2);
1726
1727 if (Min_chi < Chisq_cut) return true;
1728 else return false;
1729}
1730
1731//==============================================================================================================================
1732// Method4: loop the largest N energetic clusters on each hemisphere, then select the combination giving the smallest chisq
1733//==============================================================================================================================
1734
1735bool CgemLineFit::Loop_MaxQ() {
1736 // modified by wxn
1737 double cl_00(0), cl_01(0), cl_11(0), cl_10(0), cl_20(0), cl_21(0);
1738 int _loop(0);
1739 double maxQ_21(0), maxQ_20(0), maxQ_11(0), maxQ_10(0), maxQ_00(0);
1740 int max_21(-1), max_20(-1), max_11(-1), max_10(-1), max_00(-1);
1741 double C_00(0), C_10(0), C_11(0), C_01(0), C_20(0), C_21(0);
1742 int iC_00(0), iC_10(0), iC_11(0), iC_01(0), iC_20(0), iC_21(0);
1743 if (run < 0) Get_MCInf();
1744 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; //
1745 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; //
1746 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; //
1747 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; //
1748 vector<double> Vec_20_layer, Vec_20_phi, Vec_20_Z, Vec_20_layerid, Vec_20_v, Vec_20_flag, Vec_20_Q, Vec_20_sheetid; //
1749 vector<double> Vec_21_layer, Vec_21_phi, Vec_21_Z, Vec_21_layerid, Vec_21_v, Vec_21_flag, Vec_21_Q, Vec_21_sheetid; //
1750
1751 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
1752 RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
1753 vector<int> L_00, L_01, L_10, L_11, L_20, L_21;
1754 vector<double> Q_00, Q_01, Q_10, Q_11, Q_20, Q_21;
1755 vector<int> ZL_00, ZL_01, ZL_10, ZL_11, ZL_20, ZL_21;
1756 set<int> ZL_00_s, ZL_01_s, ZL_10_s, ZL_11_s, ZL_20_s, ZL_21_s;
1757
1758 int ic(0);
1759 if (debug) cout << "start Loop_MaxQ method" << endl;
1760
1761 for (; iter != recCgemClusterCol->end(); iter++) {
1762 _loop++;
1763 RecCgemCluster *cluster = (*iter);
1764
1765 int flag = cluster->getflag();
1766 if (flag != 2) continue;
1767 int layerid = cluster->getlayerid();
1768 double _Q = cluster->getenergydeposit();
1769 double _phi = cluster->getrecphi();
1770 double _v = cluster->getrecv();
1771 int sheetid = cluster->getsheetid();
1772 int clusterid = cluster->getclusterid();
1773 int clusterflagb = cluster->getclusterflagb();
1774 int clusterflage = cluster->getclusterflage();
1775 double x = R_layer[layerid] * cos(_phi);
1776 double y = R_layer[layerid] * sin(_phi);
1777 double z = cluster->getRecZ();
1778
1779 // --- cuts Q>10fC, size>1
1780 //RecCgemCluster* clusterX=*(recCgemClusterCol->begin()+clusterflagb);
1781 //RecCgemCluster* clusterV=*(recCgemClusterCol->begin()+clusterflage);
1782 //double QX=clusterX->getenergydeposit();
1783 //if(QX<10.) continue;
1784 //double QV=clusterV->getenergydeposit();
1785 //if(QV<10.) continue;
1786 //int size_X = 1+abs(clusterX->getclusterflage()-clusterX->getclusterflagb());
1787 //if(size_X<2) continue;
1788 //int size_V = 1+abs(clusterV->getclusterflage()-clusterV->getclusterflagb());
1789 //if(size_V<2) continue;
1790
1791 if (debug) cout << "layerid is " << layerid << endl;
1792 if (_loop > MAX_Clusters) break;
1793 ic++;
1794 if (layerid == 2 && sheetid == 1) {
1795 ZL_21_s.insert(clusterflagb);
1796 // is_cluster_exsited[0]++;
1797 }
1798 if (layerid == 2 && sheetid == 0) {
1799 ZL_20_s.insert(clusterflagb);
1800 // is_cluster_exsited[1]++;
1801 }
1802 if (layerid == 1 && sheetid == 1) {
1803 ZL_11_s.insert(clusterflagb);
1804 // is_cluster_exsited[2]++;
1805 }
1806 if (layerid == 1 && sheetid == 0) {
1807 ZL_10_s.insert(clusterflagb);
1808 // is_cluster_exsited[3]++;
1809 }
1810 if (layerid == 0 && sheetid == 0 && _phi > 0) {
1811 ZL_01_s.insert(clusterflagb);
1812 // is_cluster_exsited[4]++;
1813 }
1814 if (layerid == 0 && sheetid == 0 && _phi < 0) {
1815 ZL_00_s.insert(clusterflagb);
1816 // is_cluster_exsited[5]++;
1817 }
1818 }
1819 ZL_21.assign(ZL_21_s.begin(), ZL_21_s.end());
1820 ZL_20.assign(ZL_20_s.begin(), ZL_20_s.end());
1821 ZL_11.assign(ZL_11_s.begin(), ZL_11_s.end());
1822 ZL_10.assign(ZL_10_s.begin(), ZL_10_s.end());
1823 ZL_01.assign(ZL_01_s.begin(), ZL_01_s.end());
1824 ZL_00.assign(ZL_00_s.begin(), ZL_00_s.end());
1825 ////
1826 // if (debug) cout << "ZL_21.size() : " << ZL_21.size() << endl;
1827 // if (is_cluster_exsited[0]>0&& is_cluster_exsited[1]>0&& is_cluster_exsited[2]>0&& is_cluster_exsited[3]>0&&
1828 // is_cluster_exsited[4]>0&& is_cluster_exsited[5] > 0)
1829 //{
1830 // nEvent_with_6cluster++;
1831 // }
1832
1833 // if (is_cluster_exsited[0]>0&& is_cluster_exsited[1]>0&& is_cluster_exsited[2]>0&& is_cluster_exsited[3]> 0)
1834 //{
1835 // nEvent_with_4cluster++;
1836 // }
1837
1838 ////
1839 iter = recCgemClusterCol->begin();
1840 _loop = 0;
1841 for (; iter != recCgemClusterCol->end(); iter++) {
1842 _loop++;
1843 RecCgemCluster *cluster = (*iter);
1844 int flag = cluster->getflag();
1845 int layerid = cluster->getlayerid();
1846 int sheetid = cluster->getsheetid();
1847 double _Q = cluster->getenergydeposit();
1848 double _phi = cluster->getrecphi();
1849 double _v = cluster->getrecv();
1850 double _Z = cluster->getRecZ();
1851 int clusterid = cluster->getclusterid();
1852 int clusterflagb = cluster->getclusterflagb();
1853 int clusterflage = cluster->getclusterflage();
1854 if (flag != 0) continue;
1855 if (layerid == 2 && sheetid == 1) {
1856 for (int i_21 = 0; i_21 < ZL_21.size(); i_21++) {
1857 if (clusterid == ZL_21[i_21]) {
1858 L_21.push_back(_loop);
1859 Q_21.push_back(_Q);
1860 }
1861 }
1862 }
1863 if (layerid == 2 && sheetid == 0) {
1864 for (int i_20 = 0; i_20 < ZL_20.size(); i_20++) {
1865 if (clusterid == ZL_20[i_20]) {
1866 L_20.push_back(_loop);
1867 Q_20.push_back(_Q);
1868 }
1869 }
1870 }
1871 if (layerid == 1 && sheetid == 1) {
1872 for (int i_11 = 0; i_11 < ZL_11.size(); i_11++) {
1873 if (clusterid == ZL_11[i_11]) {
1874 L_11.push_back(_loop);
1875 Q_11.push_back(_Q);
1876 }
1877 }
1878 }
1879 if (layerid == 1 && sheetid == 0) {
1880 for (int i_10 = 0; i_10 < ZL_10.size(); i_10++) {
1881 if (clusterid == ZL_10[i_10]) {
1882 L_10.push_back(_loop);
1883 Q_10.push_back(_Q);
1884 }
1885 }
1886 }
1887 if (layerid == 0 && sheetid == 0 && _phi > 0) {
1888 for (int i_01 = 0; i_01 < ZL_01.size(); i_01++) {
1889 if (clusterid == ZL_01[i_01]) {
1890 L_01.push_back(_loop);
1891 Q_01.push_back(_Q);
1892 }
1893 }
1894 }
1895 if (layerid == 0 && sheetid == 0 && _phi < 0) {
1896 for (int i_00 = 0; i_00 < ZL_00.size(); i_00++) {
1897 if (clusterid == ZL_00[i_00]) {
1898 L_00.push_back(_loop);
1899 Q_00.push_back(_Q);
1900 }
1901 }
1902 }
1903 }
1904 // chose the N X clusters with the largest charge and save corresponding
1905 // index
1906 vector<int> L_00_s, L_01_s, L_10_s, L_11_s, L_20_s, L_21_s;
1907
1908 L_00_s = GetNMaxQ(Q_00, L_00, Nmax);
1909 L_01_s = GetNMaxQ(Q_01, L_01, Nmax);
1910 L_10_s = GetNMaxQ(Q_10, L_10, Nmax);
1911 L_11_s = GetNMaxQ(Q_11, L_11, Nmax);
1912 L_20_s = GetNMaxQ(Q_20, L_20, Nmax);
1913 L_21_s = GetNMaxQ(Q_21, L_21, Nmax);
1914 ////
1915 if (debug) cout << "L_21_s.size() : " << L_21_s.size() << endl;
1916 ////
1917
1918 if (TEST_N == 1) {
1919 L_21_s.clear();
1920 L_21_s.push_back(-1);
1921 } else if (TEST_N == 2) {
1922 L_11_s.clear();
1923 L_11_s.push_back(-1);
1924 } else if (TEST_N == 3) {
1925 L_01_s.clear();
1926 L_01_s.push_back(-1);
1927 } else if (TEST_N == 4) {
1928 L_00_s.clear();
1929 L_00_s.push_back(-1);
1930 } else if (TEST_N == 5) {
1931 L_10_s.clear();
1932 L_10_s.push_back(-1);
1933 } else if (TEST_N == 6) {
1934 L_20_s.clear();
1935 L_20_s.push_back(-1);
1936 }
1937
1938 if (debug) cout << "finish the first loop" << endl;
1939
1940 double Min_chi(1E10);
1941 NXComb = L_00_s.size() * L_01_s.size() * L_10_s.size() * L_11_s.size() * L_20_s.size() * L_21_s.size();
1942
1943 int com0, com1, com2, com3, com4, com5;
1944 for (int i_21 = 0; i_21 < L_21_s.size(); i_21++) {
1945 for (int i_20 = 0; i_20 < L_20_s.size(); i_20++) {
1946 for (int i_11 = 0; i_11 < L_11_s.size(); i_11++) {
1947 for (int i_10 = 0; i_10 < L_10_s.size(); i_10++) {
1948 for (int i_00 = 0; i_00 < L_00_s.size(); i_00++) {
1949 for (int i_01 = 0; i_01 < L_01_s.size(); i_01++) {
1950
1951 Vec_layer.clear();
1952 Vec_phi.clear();
1953 Vec_Z.clear();
1954 Vec_v.clear();
1955 Vec_layerid.clear();
1956 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
1957 Vec_flag.clear();
1958 Vec_Q.clear();
1959 Vec_sheetid.clear();
1960 Vec_m_layer.clear();
1961 Vec_m_phi.clear();
1962 Vec_m_Z.clear();
1963 Vec_m_v.clear();
1964 Vec_m_layerid.clear();
1965 Vec_m_flag.clear();
1966 Vec_m_Q.clear();
1967 Vec_m_sheetid.clear();
1968 _loop = 0;
1969 iter = recCgemClusterCol->begin();
1970
1971 vector<int> Com;
1972 Com.clear();
1973 for (; iter != recCgemClusterCol->end(); iter++) {
1974 _loop++;
1975 RecCgemCluster *cluster = (*iter);
1976 int flag = cluster->getflag();
1977 /////
1978 // if (debug) cout << "(wxn)cluster flag is " << flag << endl;
1979 /////
1980 if (!(_loop == L_21_s[i_21] || _loop == L_20_s[i_20] || _loop == L_11_s[i_11] || _loop == L_10_s[i_10] ||
1981 _loop == L_00_s[i_00] || _loop == L_01_s[i_01]))
1982 continue;
1983 /////
1984 // if (debug) cout << "(wxn)_loop is " << _loop << endl;
1985 /////
1986
1987 int layerid = cluster->getlayerid();
1988 int sheetid = cluster->getsheetid();
1989 double _phi = cluster->getrecphi();
1990 double _v = cluster->getrecv();
1991 double _Z = cluster->getRecZ();
1992 double _Q = cluster->getenergydeposit();
1993 double t_phi = cluster->getrecphi_mTPC();
1994 double t_v = cluster->getrecv_mTPC();
1995 double t_Z = cluster->getRecZ_mTPC();
1996 double merged_phi = cluster->get_merge_phi();
1997 double merged_v = cluster->get_merge_v();
1998 /////
1999 // cout << "(wxn)mergedv is " << merged_v << endl;
2000 /////
2001
2002 int vir_layerid = GetVirLay(layerid, _phi);
2003
2004 Vec_flag.push_back(flag);
2005 Vec_layerid.push_back(layerid);
2006 Vec_Virlayerid.push_back(vir_layerid);
2007 Vec_sheetid.push_back(sheetid);
2008 Vec_layer.push_back(R_layer[layerid]);
2009 Com.push_back(_loop);
2010
2011 Vec_m_phi.push_back(t_phi);
2012 Vec_m_v.push_back(t_v);
2013 Vec_m_Z.push_back(t_Z);
2014
2015 if (Switch_CCmTPC == 0) {
2016 Vec_phi.push_back(_phi);
2017 Vec_v.push_back(_v);
2018 Vec_Z.push_back(_Z);
2019 if (debug) cout << "1st phi, v, by merge: " << merged_phi << ", " << merged_v << endl;
2020 if (debug) cout << "1st phi, v, by CC: " << _phi << ", " << _v << endl;
2021 } else if (Switch_CCmTPC == 1) {
2022 Vec_phi.push_back(t_phi);
2023 Vec_v.push_back(t_v);
2024 Vec_Z.push_back(t_Z);
2025 } else if (Switch_CCmTPC == 2) {
2026
2027 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
2028 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
2029 Vec_phi.push_back(merged_phi);
2030 Vec_v.push_back(merged_v);
2031 Vec_Z.push_back(merged_Z);
2032 if (debug) cout << "1st phi, v, by merge: " << merged_phi << ", " << merged_v << endl;
2033 if (debug) cout << "1st phi, v, by CC: " << _phi << ", " << _v << endl;
2034 }
2035
2036 } // one combo
2037 if (Vec_layerid.size() != 6 && TEST_N == 0) return false;
2038 else if (Vec_layerid.size() != 5 && TEST_N > 0) return false;
2039 if (debug) cout << "before OrderCluster" << endl;
2040 OrderClusters();
2041 // 0 1 2 3
2042 // order is layer 1 , top, down, layer 0, top down
2043 // ????? How to order layer 2 ?????
2044 // OrderClusters() need modified !!!!!
2045 // 0 1 2 3 4 5
2046 // order is layer 2 , top, down, layer 1, top down, layer 0, top down
2047
2048 StraightLine *l_outer_tmp;
2049
2050 if (debug) cout << "before IniPar" << endl;
2051 // lay0top virlay
2052 // lay0down virlay
2053 if (TEST_N == 1 || TEST_N == 6) l_outer_tmp = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
2054 else if (TEST_N == 2 || TEST_N == 3 || TEST_N == 4 || TEST_N == 5)
2055 l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
2056 else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
2057
2058 if (debug) cout << "before fit" << endl;
2059 TMinuit *_fit = Fit(l_outer_tmp, fa, 1);
2060
2061 if (debug) cout << "chisq is : " << _fit->fAmin << endl;
2062 if (debug) cout << "===========================================" << endl;
2063
2064 double _chi = _fit->fAmin;
2065 if (debug) cout << "finish the fit and the chisq is " << _chi << endl;
2066 if (debug) cout << "----------------------------------------" << endl;
2067 if (_chi < Min_chi) {
2068 Min_chi = _chi;
2069 if (TEST_N == 0) C_00 = Vec_phi[5];
2070 else C_00 = -99; // Why set C_00 to -99 Guoaq--2020/11/01
2071 C_01 = Vec_phi[4]; // It is because there is only 3 clusters
2072 // selected if TEST_N!=0
2073 C_10 = Vec_phi[3];
2074 C_11 = Vec_phi[2];
2075 C_20 = Vec_phi[1];
2076 C_21 = Vec_phi[0];
2077 com0 = Com[0];
2078 com1 = Com[1];
2079 com2 = Com[2];
2080 com3 = Com[3];
2081 com4 = Com[4];
2082 com5 = Com[5];
2083 }
2084 delete _fit;
2085 delete l_outer_tmp;
2086 }
2087 }
2088 }
2089 }
2090 }
2091 }
2092
2093 if (debug) cout << "finish the second loop" << endl;
2094
2095 iter = recCgemClusterCol->begin();
2096 ZL_00.clear();
2097 ZL_01.clear();
2098 ZL_10.clear();
2099 ZL_11.clear();
2100 ZL_20.clear();
2101 ZL_21.clear();
2102 L_00.clear();
2103 L_01.clear();
2104 L_10.clear();
2105 L_11.clear();
2106 L_20.clear();
2107 L_21.clear();
2108 L_00_s.clear();
2109 L_01_s.clear();
2110 L_10_s.clear();
2111 L_11_s.clear();
2112 L_20_s.clear();
2113 L_21_s.clear();
2114 Q_00.clear();
2115 Q_01.clear();
2116 Q_10.clear();
2117 Q_11.clear();
2118 Q_20.clear();
2119 Q_21.clear();
2120 _loop = 0;
2121 int ic_tot(0);
2122
2123 if (debug) cout << "c00, c01, c11, c10 : " << C_00 << ", " << C_01 << ", " << C_11 << ", " << C_10 << endl;
2124 if (debug) cout << "ids is : " << com0 << ", " << com1 << ", " << com2 << ", " << com3 << ", " << com4 << ", " << com5 << endl;
2125 for (; iter != recCgemClusterCol->end(); iter++) {
2126 _loop++;
2127 RecCgemCluster *cluster = (*iter);
2128 int flag = cluster->getflag();
2129 if (flag != 2) continue;
2130 ic_tot++;
2131 int layerid = cluster->getlayerid();
2132 double _Q = cluster->getenergydeposit();
2133 // double t_phi=cluster->getrecphi_mTPC();
2134
2135 double _phi = 99;
2136 double _Z = cluster->getRecZ();
2137 if (Switch_CCmTPC == 0) _phi = cluster->getrecphi();
2138 else if (Switch_CCmTPC == 1) _phi = cluster->getrecphi_mTPC();
2139 else if (Switch_CCmTPC == 2) {
2140 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
2141 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
2142 _phi = cluster_x->get_merge_phi();
2143 }
2144 if (debug)
2145 cout << "2D: phi_cc, phi_merge, index " << cluster->getrecphi() << ", " << cluster->get_merge_phi() << ", " << _loop
2146 << endl;
2147 int sheetid = cluster->getsheetid();
2148 if (_Q < MinQ_Clus2D) continue; // apply a requireent on the Q of 2D cluster
2149
2150 // --- cuts Q>10fC, size>1
2151 //int clusterflage = cluster->getclusterflage();
2152 //RecCgemCluster* clusterV=*(recCgemClusterCol->begin()+clusterflage);
2153 //double QV=clusterV->getenergydeposit();
2154 //if(QV<10.) continue;
2155 //int size_V = 1+abs(clusterV->getclusterflage()-clusterV->getclusterflagb());
2156 //if(size_V<2) continue;
2157
2158 if (fabs(C_00 - _phi) < 0.001 || fabs(C_01 - _phi) < 0.001 || fabs(C_11 - _phi) < 0.001 || fabs(C_10 - _phi) < 0.001 ||
2159 fabs(C_21 - _phi) < 0.001 || fabs(C_20 - _phi) < 0.001) { // is it dangerious to compare two float values ?
2160 if (layerid == 2 && sheetid == 1) {
2161 L_21.push_back(_loop);
2162 Q_21.push_back(_Q);
2163 }
2164 if (layerid == 2 && sheetid == 0) {
2165 L_20.push_back(_loop);
2166 Q_20.push_back(_Q);
2167 }
2168 if (layerid == 1 && sheetid == 1) {
2169 L_11.push_back(_loop);
2170 Q_11.push_back(_Q);
2171 }
2172 if (layerid == 1 && sheetid == 0) {
2173 L_10.push_back(_loop);
2174 Q_10.push_back(_Q);
2175 }
2176 if (layerid == 0 && sheetid == 0 && _phi > 0) {
2177 L_01.push_back(_loop);
2178 Q_01.push_back(_Q);
2179 }
2180 if (layerid == 0 && sheetid == 0 && _phi < 0) {
2181 L_00.push_back(_loop);
2182 Q_00.push_back(_Q);
2183 }
2184 }
2185 }
2186
2187 if (debug) cout << "finish the third loop" << endl;
2188
2189 Min_chi = 1E10;
2190
2191 clst_21 = L_21.size();
2192 clst_11 = L_11.size();
2193 clst_01 = L_01.size();
2194 clst_20 = L_20.size();
2195 clst_10 = L_10.size();
2196 clst_00 = L_00.size();
2197
2198 if (debug)
2199 cout << "clst_21, clst_11, clst_01, clst_20, clst_10, clst_00 : " << clst_21 << ", " << clst_11 << ", " << clst_01 << ", "
2200 << clst_20 << ", " << clst_10 << ", " << clst_00 << endl;
2201
2202 L_00_s = GetNMaxQ(Q_00, L_00, Nmax);
2203 L_01_s = GetNMaxQ(Q_01, L_01, Nmax);
2204 L_10_s = GetNMaxQ(Q_10, L_10, Nmax);
2205 L_11_s = GetNMaxQ(Q_11, L_11, Nmax);
2206 L_20_s = GetNMaxQ(Q_20, L_20, Nmax);
2207 L_21_s = GetNMaxQ(Q_21, L_21, Nmax);
2208
2209 if (TEST_N == 1) {
2210 L_21_s.clear();
2211 L_21_s.push_back(-1);
2212 } else if (TEST_N == 2) {
2213 L_11_s.clear();
2214 L_11_s.push_back(-1);
2215 } else if (TEST_N == 3) {
2216 L_01_s.clear();
2217 L_01_s.push_back(-1);
2218 } else if (TEST_N == 4) {
2219 L_00_s.clear();
2220 L_00_s.push_back(-1);
2221 } else if (TEST_N == 5) {
2222 L_10_s.clear();
2223 L_10_s.push_back(-1);
2224 } else if (TEST_N == 6) {
2225 L_20_s.clear();
2226 L_20_s.push_back(-1);
2227 }
2228
2229 NVComb = L_00_s.size() * L_01_s.size() * L_10_s.size() * L_11_s.size() * L_20_s.size() * L_21_s.size();
2230
2231 for (int i_21 = 0; i_21 < L_21_s.size(); i_21++) {
2232 for (int i_20 = 0; i_20 < L_20_s.size(); i_20++) {
2233 for (int i_11 = 0; i_11 < L_11_s.size(); i_11++) {
2234 for (int i_10 = 0; i_10 < L_10_s.size(); i_10++) {
2235 for (int i_00 = 0; i_00 < L_00_s.size(); i_00++) {
2236 for (int i_01 = 0; i_01 < L_01_s.size(); i_01++) {
2237
2238 Vec_layer.clear();
2239 Vec_phi.clear();
2240 Vec_Z.clear();
2241 Vec_v.clear();
2242 Vec_layerid.clear();
2243 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
2244 Vec_flag.clear();
2245 Vec_Q.clear();
2246 Vec_sheetid.clear();
2247 Vec_m_layer.clear();
2248 Vec_m_phi.clear();
2249 Vec_m_Z.clear();
2250 Vec_m_v.clear();
2251 Vec_m_layerid.clear();
2252 Vec_m_flag.clear();
2253 Vec_m_Q.clear();
2254 Vec_m_sheetid.clear();
2255 _loop = 0;
2256 iter = recCgemClusterCol->begin();
2257 for (; iter != recCgemClusterCol->end(); iter++) {
2258 _loop++;
2259 RecCgemCluster *cluster = (*iter);
2260 int flag = cluster->getflag();
2261 if (flag != 2) continue;
2262 if (!(_loop == L_21_s[i_21] || _loop == L_20_s[i_20] || _loop == L_11_s[i_11] || _loop == L_10_s[i_10] ||
2263 _loop == L_00_s[i_00] || _loop == L_01_s[i_01]))
2264 continue;
2265
2266 int layerid = cluster->getlayerid();
2267 int sheetid = cluster->getsheetid();
2268 double _phi = cluster->getrecphi();
2269 double _v = cluster->getrecv();
2270 double _Z = cluster->getRecZ();
2271 double _Q = cluster->getenergydeposit();
2272 double t_phi = cluster->getrecphi_mTPC();
2273 double t_v = cluster->getrecv_mTPC();
2274 double t_Z = cluster->getRecZ_mTPC();
2275
2276 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
2277 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
2278 double merged_phi = cluster_x->get_merge_phi();
2279 RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
2280 double merged_v = cluster_v->get_merge_v();
2281
2282 int vir_layerid = GetVirLay(layerid, _phi);
2283
2284 Vec_flag.push_back(flag);
2285 Vec_layerid.push_back(layerid);
2286 Vec_sheetid.push_back(sheetid);
2287 Vec_layer.push_back(R_layer[layerid]);
2288 Vec_Virlayerid.push_back(vir_layerid);
2289
2290 Vec_m_phi.push_back(t_phi);
2291 Vec_m_v.push_back(t_v);
2292 Vec_m_Z.push_back(t_Z);
2293
2294 if (Switch_CCmTPC == 0) {
2295 Vec_phi.push_back(_phi);
2296 Vec_v.push_back(_v);
2297 Vec_Z.push_back(_Z);
2298 } else if (Switch_CCmTPC == 1) {
2299 Vec_phi.push_back(t_phi);
2300 Vec_v.push_back(t_v);
2301 Vec_Z.push_back(t_Z);
2302 } else if (Switch_CCmTPC == 2) {
2303
2304 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
2305 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
2306 Vec_phi.push_back(merged_phi);
2307 Vec_v.push_back(merged_v);
2308 Vec_Z.push_back(merged_Z);
2309 }
2310
2311 } // one combo
2312 if (Vec_layerid.size() != 6 && TEST_N == 0) return false;
2313 else if (Vec_layerid.size() != 5 && TEST_N > 0) return false;
2314 OrderClusters();
2315
2316 // discard the combination if V strip is re-used !
2317 if (Vec_v[0] == Vec_v[1] || Vec_v[0] == Vec_v[2] || Vec_v[0] == Vec_v[3] || Vec_v[0] == Vec_v[4] ||
2318 Vec_v[0] == Vec_v[5] || Vec_v[1] == Vec_v[2] || Vec_v[1] == Vec_v[3] || Vec_v[1] == Vec_v[4] ||
2319 Vec_v[1] == Vec_v[5] || Vec_v[2] == Vec_v[3] || Vec_v[2] == Vec_v[4] || Vec_v[2] == Vec_v[5] ||
2320 Vec_v[3] == Vec_v[4] || Vec_v[3] == Vec_v[5] || Vec_v[4] == Vec_v[5])
2321 continue;
2322 StraightLine *l_outer_tmp;
2323
2324 if (TEST_N == 1 || TEST_N == 6) l_outer_tmp = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
2325 else if (TEST_N == 2 || TEST_N == 3 || TEST_N == 4 || TEST_N == 5)
2326 l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
2327 else if (TEST_N == 0) l_outer_tmp = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
2328
2329 TMinuit *_fit = Fit(l_outer_tmp, fa, 2);
2330 double _chi = _fit->fAmin;
2331 if (_chi < Min_chi) {
2332 Min_chi = _chi;
2333 iC_00 = L_00_s[i_00];
2334 iC_01 = L_01_s[i_01];
2335 iC_10 = L_10_s[i_10];
2336 iC_11 = L_11_s[i_11];
2337 iC_20 = L_20_s[i_20];
2338 iC_21 = L_21_s[i_21];
2339 }
2340 delete _fit;
2341 delete l_outer_tmp;
2342 }
2343 }
2344 }
2345 }
2346 }
2347 }
2348
2349 Vec_layer.clear();
2350 Vec_phi.clear();
2351 Vec_Z.clear();
2352 Vec_v.clear();
2353 Vec_layerid.clear();
2354 Vec_Virlayerid.clear(); // Add by Guoaq -- 2020/11/01
2355 Vec_flag.clear();
2356 Vec_Q.clear();
2357 Vec_sheetid.clear();
2358 Vec_xCluSize.clear();
2359 Vec_vCluSize.clear();
2360 Vec_XQ.clear();
2361 Vec_VQ.clear();
2362 Vec_XstripQ.clear();
2363 Vec_VstripQ.clear();
2364 Vec_XstripT.clear();
2365 Vec_VstripT.clear();
2366 Vec_XstripTf.clear();
2367 Vec_VstripTf.clear();
2368 Vec_XstripID.clear();
2369 Vec_VstripID.clear();
2370 Vec_m_layer.clear();
2371 Vec_m_phi.clear();
2372 Vec_m_Z.clear();
2373 Vec_m_v.clear();
2374 Vec_m_layerid.clear();
2375 Vec_m_flag.clear();
2376 Vec_m_Q.clear();
2377 Vec_m_sheetid.clear();
2378 if (debug) cout << "finish the forth loop" << endl;
2379 _loop = 0;
2380 ic = 0;
2381 iter = recCgemClusterCol->begin();
2382 for (; iter != recCgemClusterCol->end(); iter++) {
2383 _loop++;
2384 RecCgemCluster *cluster = (*iter);
2385 int flag = cluster->getflag();
2386 int layerid = cluster->getlayerid();
2387 int sheetid = cluster->getsheetid();
2388 double _phi = cluster->getrecphi();
2389 double _v = cluster->getrecv();
2390 double _Z = cluster->getRecZ();
2391 double _Q = cluster->getenergydeposit();
2392 double t_phi = cluster->getrecphi_mTPC();
2393 double t_v = cluster->getrecv_mTPC();
2394 double t_Z = cluster->getRecZ_mTPC();
2395
2396 int clusterid = cluster->getclusterid();
2397 // int clusterflagb=cluster->getclusterflagb();
2398 // int clusterflage=cluster->getclusterflage();
2399 double x = R_layer[layerid] * cos(_phi);
2400 double y = R_layer[layerid] * sin(_phi);
2401 double z = cluster->getRecZ();
2402 int vir_layerid = GetVirLay(layerid, _phi);
2403 if (_loop > MAX_Clusters) break;
2404 if (flag != 2) continue;
2405
2406 RecCgemClusterCol::iterator iter_begin = recCgemClusterCol->begin();
2407 RecCgemCluster *cluster_x = *(iter_begin + cluster->getclusterflagb());
2408 double merged_phi = cluster_x->get_merge_phi();
2409 RecCgemCluster *cluster_v = *(iter_begin + cluster->getclusterflage());
2410 double merged_v = cluster_v->get_merge_v();
2411
2412 tr_id_cluster[ic] = clusterid;
2413 tr_x_cluster[ic] = x;
2414 tr_y_cluster[ic] = y;
2415 tr_z_cluster[ic] = z;
2416 tr_Q_cluster[ic] = _Q;
2417 tr_phi_cluster[ic] = _phi;
2418 tr_layer_cluster[ic] = layerid;
2419 tr_sheet_cluster[ic] = sheetid;
2420 tr_Is_selected_cluster[ic] = 0;
2421 // if ((_loop == iC_00 || _loop == iC_10 || _loop == iC_20 || _loop == iC_01 || _loop == iC_11 || _loop == iC_21))
2422 // {tr_Is_selected_cluster[ic] = 1;}
2423
2424 // cout << "Test : ic = " << ic << " tr_layer_cluster = " << tr_layer_cluster[ic] << " tr_sheet_cluster = " <<
2425 // tr_sheet_cluster[ic] << " tr_phi_cluster = " << tr_phi_cluster[ic] << " tr_Is_selected_cluster = " <<
2426 // tr_Is_selected_cluster[ic] << endl;
2427 ic++;
2428 if (!(_loop == iC_00 || _loop == iC_10 || _loop == iC_20 || _loop == iC_01 || _loop == iC_11 || _loop == iC_21)) continue;
2429 tr_Is_selected_cluster[ic - 1] = 1;
2430 SmartDataPtr<CgemDigiCol> cgemDigiCol(eventSvc(), "/Event/Digi/CgemDigiCol");
2431 int clusterflagb = cluster->getclusterflagb();
2432 int clusterflage = cluster->getclusterflage();
2433 RecCgemClusterCol::iterator id_XCluster = recCgemClusterCol->begin() + clusterflagb;
2434 RecCgemCluster *xcluster = (*id_XCluster);
2435 int xclustersize = abs(xcluster->getclusterflage() - xcluster->getclusterflagb()) + 1;
2436 double X_Q = xcluster->getenergydeposit();
2437 RecCgemClusterCol::iterator id_VCluster = recCgemClusterCol->begin() + clusterflage;
2438 RecCgemCluster *vcluster = (*id_VCluster);
2439 int vclustersize = abs(vcluster->getclusterflage() - vcluster->getclusterflagb()) + 1;
2440 double V_Q = vcluster->getenergydeposit();
2441 int xclusterb = xcluster->getclusterflagb();
2442 int xclustere = xcluster->getclusterflage();
2443 int vclusterb = vcluster->getclusterflagb();
2444 int vclustere = vcluster->getclusterflage();
2445 int layerid_str, sheetid_str, stripid_str, time_chnnel;
2446 int nxstrips = 0;
2447 int nvstrips = 0;
2448 double energydeposit;
2449 double Q_fC, T_ns, Tf_ns;
2450 CgemDigiCol::iterator iter_digi = cgemDigiCol->begin();
2451 Vec_x_stripsQ.clear();
2452 Vec_v_stripsQ.clear();
2453 Vec_x_stripsT.clear();
2454 Vec_v_stripsT.clear();
2455 Vec_x_stripsTf.clear();
2456 Vec_v_stripsTf.clear();
2457 Vec_x_stripsid.clear();
2458 Vec_v_stripsid.clear();
2459 for (; iter_digi != cgemDigiCol->end(); iter_digi++) {
2460 layerid_str = CgemID::layer((*iter_digi)->identify());
2461 sheetid_str = CgemID::sheet((*iter_digi)->identify());
2462 stripid_str = CgemID::strip((*iter_digi)->identify());
2463 bool striptype = CgemID::is_xstrip((*iter_digi)->identify());
2464 int flagxv;
2465 if (striptype == true) flagxv = 0;
2466 else flagxv = 1;
2467 if (layerid == layerid_str && sheetid == sheetid_str) {
2468 Q_fC = (*iter_digi)->getCharge_fc();
2469 T_ns = (*iter_digi)->getTime_ns();
2470 Tf_ns = get_Time(iter_digi);
2471 if (flagxv == 0) {
2472 if (stripid_str >= xclusterb && stripid_str <= xclustere) {
2473 // cout<<"stripid_str = "<<stripid_str<<" xclusterb =
2474 // "<<xclusterb<<" xclustere = "<<xclustere<<" "<<flagxv<<endl;
2475 if (xclustersize < 100) {
2476 Vec_x_stripsid.push_back(stripid_str);
2477 Vec_x_stripsQ.push_back(Q_fC);
2478 Vec_x_stripsT.push_back(T_ns);
2479 Vec_x_stripsTf.push_back(Tf_ns);
2480 nxstrips++;
2481 }
2482 }
2483 }
2484 if (flagxv == 1) {
2485 if (stripid_str >= vclusterb && stripid_str <= vclustere) {
2486 if (vclustersize < 100) {
2487 Vec_v_stripsid.push_back(stripid_str);
2488 Vec_v_stripsQ.push_back(Q_fC);
2489 Vec_v_stripsT.push_back(T_ns);
2490 Vec_v_stripsTf.push_back(Tf_ns);
2491 nvstrips++;
2492 }
2493 }
2494 }
2495 }
2496 }
2497
2498 Vec_flag.push_back(flag);
2499 Vec_layerid.push_back(layerid);
2500 Vec_Virlayerid.push_back(vir_layerid);
2501 Vec_sheetid.push_back(sheetid);
2502 Vec_layer.push_back(R_layer[layerid]);
2503 Vec_xCluSize.push_back(xclustersize);
2504 Vec_vCluSize.push_back(vclustersize);
2505 Vec_XQ.push_back(X_Q);
2506 Vec_VQ.push_back(V_Q);
2507 Vec_XstripQ.push_back(Vec_x_stripsQ);
2508 Vec_VstripQ.push_back(Vec_v_stripsQ);
2509 Vec_XstripT.push_back(Vec_x_stripsT);
2510 Vec_VstripT.push_back(Vec_v_stripsT);
2511 Vec_XstripTf.push_back(Vec_x_stripsTf);
2512 Vec_VstripTf.push_back(Vec_v_stripsTf);
2513 Vec_XstripID.push_back(Vec_x_stripsid);
2514 Vec_VstripID.push_back(Vec_v_stripsid);
2515
2516 Vec_m_phi.push_back(t_phi);
2517 Vec_m_v.push_back(t_v);
2518 Vec_m_Z.push_back(t_Z);
2519
2520 if (Switch_CCmTPC == 0) {
2521 Vec_phi.push_back(_phi);
2522 Vec_v.push_back(_v);
2523 Vec_Z.push_back(_Z);
2524 } else if (Switch_CCmTPC == 1) {
2525 Vec_phi.push_back(t_phi);
2526 Vec_v.push_back(t_v);
2527 Vec_Z.push_back(t_Z);
2528 } else if (Switch_CCmTPC == 2) {
2529
2530 CgemGeoReadoutPlane *A_readoutPlane = m_cgemGeomSvc->getReadoutPlane(layerid, sheetid);
2531 double merged_Z = A_readoutPlane->getZFromPhiV(merged_phi, merged_v);
2532 Vec_phi.push_back(merged_phi);
2533 Vec_v.push_back(merged_v);
2534 Vec_Z.push_back(merged_Z);
2535 }
2536
2537 cluster->setTrkId(0);
2538 vecclusters.push_back(cluster);
2539
2540 } // one combo
2541 NC = vecclusters.size();
2542 tr_ncluster = ic;
2543 if (debug) cout << "finish the fifth loop" << endl;
2544
2545 for (int i_cl = 0; i_cl < Vec_layerid.size(); i_cl++) {
2546 if (Vec_layerid[i_cl] == 2) _clst_2++;
2547 if (Vec_layerid[i_cl] == 1) _clst_1++;
2548 if (Vec_layerid[i_cl] == 0) _clst_0++;
2549 }
2550 if (debug) cout << "Vec_layerid.size is " << Vec_layerid.size() << endl;
2551
2552 if (Vec_layerid.size() != 6 && TEST_N == 0) return false;
2553 else if (Vec_layerid.size() != 5 && TEST_N > 0) return false;
2554 // OrderClusters();
2555 OrderClusterSizeQ();
2556
2557 if (TEST_N == 1 || TEST_N == 6) l_outer = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
2558 else if (TEST_N == 2 || TEST_N == 3 || TEST_N == 4 || TEST_N == 5)
2559 l_outer = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
2560 else if (TEST_N == 0) l_outer = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
2561 if (debug) {
2562 cout << "Event " << event << ", Clu. ID are " << iC_00 << ", " << iC_01 << ", " << iC_10 << ", " << iC_11 << ", " << iC_20
2563 << ", " << iC_21 << ", Loop_MaxQ Min_chi = " << Min_chi << endl;
2564 }
2565
2566 if (Min_chi < Chisq_cut) return true;
2567 else return false;
2568}
2569
2570//========================================================================================================
2571// useful functions
2572//========================================================================================================
2573
2574StatusCode CgemLineFit ::registerTrack(RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol) {
2575 MsgStream log(msgSvc(), name());
2576 StatusCode sc;
2577 IDataProviderSvc *eventSvc = NULL;
2578 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
2579 if (eventSvc) {
2580 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
2581 } else {
2582 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
2583 return StatusCode::FAILURE;
2584 // return StatusCode::SUCCESS;
2585 }
2586 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc *>(eventSvc);
2587 ;
2588
2589 // --- register RecMdcTrack
2590 recMdcTrackCol = new RecMdcTrackCol;
2591 DataObject *aRecMdcTrackCol;
2592 eventSvc->findObject("/Event/Recon/RecMdcTrackCol", aRecMdcTrackCol);
2593 if (aRecMdcTrackCol != NULL) {
2594 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
2595 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
2596 }
2597
2598 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
2599 if (sc.isFailure()) {
2600 log << MSG::FATAL << " Could not register RecMdcTrack collection" << endreq;
2601 return StatusCode::FAILURE;
2602 }
2603 log << MSG::INFO << "RecMdcTrackCol registered successfully!" << endreq;
2604
2605 recMdcHitCol = new RecMdcHitCol;
2606 DataObject *aRecMdcHitCol;
2607 eventSvc->findObject("/Event/Recon/RecMdcHitCol", aRecMdcHitCol);
2608 if (aRecMdcHitCol != NULL) {
2609 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
2610 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
2611 }
2612
2613 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
2614 if (sc.isFailure()) {
2615 log << MSG::FATAL << " Could not register RecMdcHit collection" << endreq;
2616 return StatusCode::FAILURE;
2617 }
2618 log << MSG::INFO << "RecMdcHitCol registered successfully!" << endreq;
2619
2620 return StatusCode::SUCCESS;
2621} // end of stroeTracks
2622
2623void CgemLineFit::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
2624 double chisq = 0;
2625 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
2626
2627 if (Vec_Virlayerid[i] == 5 && bool(sheet_flag & f21)) continue;
2628 if (Vec_Virlayerid[i] == 4 && bool(sheet_flag & f20)) continue;
2629 if (Vec_Virlayerid[i] == 3 && bool(sheet_flag & f11)) continue;
2630 if (Vec_Virlayerid[i] == 2 && bool(sheet_flag & f10)) continue;
2631 if (Vec_Virlayerid[i] == 1 && bool(sheet_flag & f01)) continue;
2632 if (Vec_Virlayerid[i] == 0 && bool(sheet_flag & f00)) continue;
2633 double dx1 = dx(Vec_Virlayerid[i], Vec_layer[i], par[0], par[1], dz_set, tanl_set, Vec_phi[i]);
2634 StraightLine *l = new StraightLine(par[0], par[1], dz_set, tanl_set);
2635 double flagg = Vec_flag[i];
2636 if (flagg == 2) flagg = 0;
2637 double _sigma2 = CalSigma2(Vec_layerid[i], flagg, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
2638 chisq += dx1 / _sigma2;
2639 delete l;
2640 }
2641 f = chisq;
2642}
2643
2644void CgemLineFit::fcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
2645 double chisq = 0;
2646
2647 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
2648 if (Vec_Virlayerid[i] == 5 && bool(sheet_flag & f21)) continue;
2649 if (Vec_Virlayerid[i] == 4 && bool(sheet_flag & f20)) continue;
2650 if (Vec_Virlayerid[i] == 3 && bool(sheet_flag & f11)) continue;
2651 if (Vec_Virlayerid[i] == 2 && bool(sheet_flag & f10)) continue;
2652 if (Vec_Virlayerid[i] == 1 && bool(sheet_flag & f01)) continue;
2653 if (Vec_Virlayerid[i] == 0 && bool(sheet_flag & f00)) continue;
2654
2655 if (Vec_flag[i] != 2) continue;
2656 double flagg = Vec_flag[i];
2657 if (flagg == 2) flagg = 1;
2658 double dv1 = dV(Vec_Virlayerid[i], Vec_layer[i], dr_set, phi0_set, par[0], par[1], Vec_phi[i], Vec_v[i]);
2659 StraightLine *l = new StraightLine(dr_set, phi0_set, par[0], par[1]);
2660 // double _sigma2 = CalSigma2(Vec_layerid[i], Vec_flag[i], Vec_phi[i], Vec_v[i], Vec_Z[i], l);
2661 double _sigma2 = CalSigma2(Vec_layerid[i], flagg, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
2662 chisq += dv1 / _sigma2;
2663 delete l;
2664 }
2665 f = chisq;
2666}
2667
2668void CgemLineFit::fcn3(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
2669
2670 double chisq = 0;
2671 double dv1(0), dx1(0);
2672 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
2673 if (Vec_flag[i] == 2) {
2674
2675 if (Vec_Virlayerid[i] == 5 && bool(sheet_flag & f21)) continue;
2676 if (Vec_Virlayerid[i] == 4 && bool(sheet_flag & f20)) continue;
2677 if (Vec_Virlayerid[i] == 3 && bool(sheet_flag & f11)) continue;
2678 if (Vec_Virlayerid[i] == 2 && bool(sheet_flag & f10)) continue;
2679 if (Vec_Virlayerid[i] == 1 && bool(sheet_flag & f01)) continue;
2680 if (Vec_Virlayerid[i] == 0 && bool(sheet_flag & f00)) continue;
2681
2682 dv1 = dV(Vec_Virlayerid[i], Vec_layer[i], par[0], par[1], par[2], par[3], Vec_phi[i], Vec_v[i]);
2683 StraightLine *l = new StraightLine(par[0], par[1], par[2], par[3]);
2684 double _sigma2 = CalSigma2(Vec_layerid[i], 1, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
2685 chisq += (dv1 / _sigma2);
2686 _sigma2 = CalSigma2(Vec_layerid[i], 0, Vec_phi[i], Vec_v[i], Vec_Z[i], l);
2687 dx1 = dx(Vec_Virlayerid[i], Vec_layer[i], par[0], par[1], par[2], par[3], Vec_phi[i]);
2688 chisq += dx1 / _sigma2;
2689 delete l;
2690 }
2691 }
2692 f = chisq;
2693}
2694
2695double CgemLineFit::dx(double vir_layerid, double R, double dr, double phi0, double dz, double tanl, double x) {
2696 if (dr > R) return 9999999999999;
2697 StraightLine l1(dr, phi0, dz, tanl);
2698 double phiV_up[2], phiV_down[2];
2699 HepPoint3D Up, Down;
2700 double phiV_up1_old[2], phiV_down1_old[2];
2701 HepPoint3D Up1_old, Down1_old;
2702 bool findinter = false;
2703 if (!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
2704 else findinter = Mp->getPointAligned(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
2705
2706 if (!findinter) return 9999999999999;
2707
2708 phiV_up[0] = fmod((phiV_up[0]), (2.0 * TMath::Pi()));
2709 phiV_down[0] = fmod((phiV_down[0]), (2.0 * TMath::Pi()));
2710
2711 if (phiV_up[0] < 0) phiV_up[0] += 2.0 * TMath::Pi();
2712 if (phiV_down[0] < 0) phiV_down[0] += 2.0 * TMath::Pi();
2713
2714 double dx1 = Min_diff(x, phiV_down[0]);
2715 double dx2 = Min_diff(x, phiV_up[0]);
2716 double mdv = min(pow(dx1, 2.0), pow(dx2, 2.0));
2717 return mdv * R * R;
2718}
2719
2720double CgemLineFit::dV(int vir_layerid, double R, double dr, double phi0, double z0, double tglam, double x, double V) {
2721
2722 HepPoint3D Up, Down;
2723 StraightLine l1(dr, phi0, z0, tglam);
2724 double phiV_up[2], phiV_down[2];
2725 bool findinter = false;
2726 if (!Align_FLAG) findinter = Mp->getPointIdealGeom(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
2727 else findinter = Mp->getPointAligned(vir_layerid, l1, Up, Down, phiV_up, phiV_down);
2728
2729 if (!findinter) return 9999999999999;
2730
2731 double dx1 = Min_diff(x, phiV_down[0]);
2732 double dx2 = Min_diff(x, phiV_up[0]);
2733
2734 if (dx1 < dx2) return pow(fabs(phiV_down[1] - V), 2.0);
2735 else return pow(fabs(phiV_up[1] - V), 2.0);
2736}
2737
2738TMinuit *CgemLineFit::Fit(StraightLine *l_outer, int _sheet_flag, int typ) {
2739 // typ == 0 : 2D fit on XY plane and SZ plane indipendently, then use the fit results as initial values to perform 3D fit!
2740 // typ == 1 : 2D fit on XY plane
2741 // typ == 2 : 3D fit on XY and SZ plane
2742
2743 double trkpar[4] = {-999, -999, -999, -999};
2744 double trkparErr[4] = {-999, -999, -999, -999};
2745
2746 sheet_flag = ~_sheet_flag;
2747 dz_set = l_outer->dz();
2748 tanl_set = l_outer->tanl();
2749 Int_t ierflg, npari, nparx, istat1, istat2, istat3;
2750 Double_t fmin, edm, errdef;
2751 Double_t arglist[10];
2752
2753 arglist[0] = 2000000;
2754 arglist[1] = 0.0001;
2755 if (typ == 2) {
2756 TMinuit *Min3 = new TMinuit(4);
2757 Min3->SetFCN(fcn3);
2758 Min3->SetPrintLevel(-1);
2759 Min3->mnparm(0, "dr", l_outer->dr(), 0.1, -1 * R_layer[2], R_layer[2], ierflg);
2760 Min3->mnparm(1, "phi0", l_outer->phi0(), 0.01, -100 * acos(-1), 100 * acos(-1),
2761 ierflg); // Why the range of phi0 is so big -- Guoaq--2020/11/01
2762 Min3->mnparm(2, "dz", l_outer->dz(), 0.1, -0.5 * length, 0.5 * length, ierflg);
2763 Min3->mnparm(3, "tanL", l_outer->tanl(), 0.01, -9999, 9999, ierflg);
2764 Min3->mnexcm("MIGRAD", arglist, 2, ierflg);
2765
2766 return Min3;
2767 }
2768
2769 if (typ == 0) {
2770
2771 arglist[0] = 2000000;
2772 arglist[1] = 0.0001;
2773 }
2774 TMinuit *Min = new TMinuit(2);
2775 Min->SetPrintLevel(-1);
2776 Min->SetFCN(fcn);
2777 Min->SetErrorDef(1.0); // For Chisq
2778 Min->mnparm(0, "dr", l_outer->dr(), 0.001, -1 * R_layer[2], R_layer[2], ierflg);
2779 Min->mnparm(1, "phi0", l_outer->phi0(), 0.001, -100 * acos(-1), 100 * acos(-1), ierflg);
2780 arglist[0] = 2000000;
2781 arglist[1] = 0.001;
2782 Min->mnexcm("MIGRAD", arglist, 2, ierflg);
2783 Min->mnstat(fmin, edm, errdef, npari, nparx, istat1);
2784 Min->GetParameter(0, trkpar[0], trkparErr[0]);
2785 // cout <<"typ1 : trkparErr[0] ----> " <<trkparErr[0] << endl;
2786 Min->GetParameter(1, trkpar[1], trkparErr[1]);
2787
2788 to_0_2pi(trkpar[1]);
2789 dr_set = trkpar[0];
2790 phi0_set = trkpar[1];
2791
2792 // cout << "trkpar[0] = " << trkpar[0] << " trkpar[1] = " << trkpar[1] << endl;
2793 chisq_1 = Min->fAmin;
2794 // cout << "chisq_1 : " << chisq_1 << endl;
2795 if (typ == 1) return Min;
2796 else delete Min;
2797
2798 /////////////////////////Fit to z0 & tanlam////////////
2799 // arglist[0] = 2000000;
2800 // arglist[1] = 0.0001;
2801 TMinuit *Min2 = new TMinuit(2);
2802 Min2->SetFCN(fcn2);
2803 Min2->SetPrintLevel(-1);
2804 Min2->mnparm(0, "dz", l_outer->dz(), 0.001, -0.5 * length, 0.5 * length, ierflg);
2805 Min2->mnparm(1, "tanL", l_outer->tanl(), 0.001, -9999, 9999, ierflg);
2806 Min2->mnexcm("MIGRAD", arglist, 2, ierflg);
2807 Min2->mnstat(fmin, edm, errdef, npari, nparx, istat2);
2808 for (int i = 0; i < 2; i++) { Min2->GetParameter(i, trkpar[i + 2], trkparErr[i + 2]); }
2809 chisq_2 = Min2->fAmin;
2810
2811 delete Min2;
2812 //////////////////re-Fit 4 parameters to get correct EDM /////////////
2813
2814 // arglist[0] = 2000000;
2815 // arglist[1] = 0.0001;
2816 TMinuit *Min3 = new TMinuit(4);
2817 Min3->SetFCN(fcn3);
2818 Min3->SetPrintLevel(-1);
2819 Min3->mnparm(0, "dr", trkpar[0], 0.001, -1 * R_layer[2], R_layer[2], ierflg);
2820 Min3->mnparm(1, "phi0", trkpar[1], 0.001, -100 * acos(-1), 100 * acos(-1), ierflg);
2821 Min3->mnparm(2, "dz", trkpar[2], 0.001, -0.5 * length, 0.5 * length, ierflg);
2822 Min3->mnparm(3, "tanL", trkpar[3], 0.001, -9999, 9999, ierflg);
2823 Min3->mnexcm("MIGRAD", arglist, 2, ierflg);
2824 /////
2825 // Min3->mnstat(fmin, edm, errdef, npari, nparx, istat3);
2826 // for (int i = 0; i < 4; i++) { Min3->GetParameter(i, trkpar[i], trkparErr[i]); }
2827 // cout <<"typ0 : trkparErr[0] ----> " <<trkparErr[0] << endl;
2828 /////
2829
2830 if (typ == 0) return Min3;
2831 else delete Min3;
2832}
2833
2834void CgemLineFit::OrderClusters() {
2835 int noswap(0);
2836 while (noswap != Vec_flag.size() - 1) {
2837 noswap = 0;
2838 for (int i = 0; i < Vec_flag.size() - 1; i++) {
2839 //cout << "In OrderClusters() i = " << i << " Vec_flag.size() = " << Vec_flag.size() << endl;
2840 if (Vec_flag[i] < Vec_flag[i + 1]) {
2841 Swap(i, i + 1);
2842 //cout << "layerid = " << Vec_layerid[i] << " sheetid = " << Vec_sheetid[i] << " flag = " << Vec_flag[i] << " phi = "
2843 //<< Vec_phi[i] << endl;
2844 } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] < Vec_layerid[i + 1]) {
2845 Swap(i, i + 1);
2846 //cout << "layerid = " << Vec_layerid[i] << " sheetid = " << Vec_sheetid[i] << " flag = " << Vec_flag[i] << " phi = "
2847 //<< Vec_phi[i] << endl;
2848 } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] == Vec_layerid[i + 1] && (Vec_phi[i] < Vec_phi[i + 1])) {
2849 Swap(i, i + 1);
2850 //cout << "layerid = " << Vec_layerid[i] << " sheetid = " << Vec_sheetid[i] << " flag = " << Vec_flag[i] << " phi = "
2851 //<< Vec_phi[i] << endl;
2852 } else {
2853 //cout << "xxxxxxx" << endl;
2854 //cout << "(i)layerid = " << Vec_layerid[i] << " sheetid = " << Vec_sheetid[i] << "Virlayerid = " << Vec_Virlayerid[i]
2855 //<< " flag = " << Vec_flag[i] << " phi = " << Vec_phi[i] << endl; cout << "(i+1)layerid = " << Vec_layerid[i+1] << "sheetid = " << Vec_sheetid[i+1] << "Virlayerid = " << Vec_Virlayerid[i+1] << " flag = " << Vec_flag[i+1] << " phi = "
2856 //<< Vec_phi[i+1] << endl;
2857 noswap++;
2858 }
2859 }
2860 }
2861}
2862
2863void CgemLineFit::Swap(int i, int j) {
2864
2865 swap(Vec_layer[i], Vec_layer[j]);
2866 swap(Vec_sheetid[i], Vec_sheetid[j]);
2867 swap(Vec_phi[i], Vec_phi[j]);
2868 swap(Vec_Z[i], Vec_Z[j]);
2869 swap(Vec_v[i], Vec_v[j]);
2870 swap(Vec_layerid[i], Vec_layerid[j]);
2871 swap(Vec_Virlayerid[i], Vec_Virlayerid[j]);
2872 swap(Vec_flag[i], Vec_flag[j]);
2873 if (Vec_m_phi.size() == Vec_flag.size() && (Method == 1 || Method == 2)) {
2874 swap(Vec_m_phi[i], Vec_m_phi[j]);
2875 swap(Vec_m_Z[i], Vec_m_Z[j]);
2876 swap(Vec_m_v[i], Vec_m_v[j]);
2877 }
2878}
2879
2880void CgemLineFit::OrderClusterSizeQ() {
2881 int noswap(0);
2882 while (noswap != Vec_flag.size() - 1) {
2883 noswap = 0;
2884 for (int i = 0; i < Vec_flag.size() - 1; i++) {
2885 if (Vec_flag[i] < Vec_flag[i + 1]) {
2886 Swap(i, i + 1);
2887 swap(Vec_xCluSize[i], Vec_xCluSize[i + 1]);
2888 swap(Vec_vCluSize[i], Vec_vCluSize[i + 1]);
2889 swap(Vec_XQ[i], Vec_XQ[i + 1]);
2890 swap(Vec_VQ[i], Vec_VQ[i + 1]);
2891 swap(Vec_XstripQ[i], Vec_XstripQ[i + 1]);
2892 swap(Vec_VstripQ[i], Vec_VstripQ[i + 1]);
2893 swap(Vec_XstripT[i], Vec_XstripT[i + 1]);
2894 swap(Vec_VstripT[i], Vec_VstripT[i + 1]);
2895 swap(Vec_XstripTf[i], Vec_XstripTf[i + 1]);
2896 swap(Vec_VstripTf[i], Vec_VstripTf[i + 1]);
2897 swap(Vec_XstripID[i], Vec_XstripID[i + 1]);
2898 swap(Vec_VstripID[i], Vec_VstripID[i + 1]);
2899 } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] < Vec_layerid[i + 1]) {
2900 Swap(i, i + 1);
2901 swap(Vec_xCluSize[i], Vec_xCluSize[i + 1]);
2902 swap(Vec_vCluSize[i], Vec_vCluSize[i + 1]);
2903 swap(Vec_XQ[i], Vec_XQ[i + 1]);
2904 swap(Vec_VQ[i], Vec_VQ[i + 1]);
2905 swap(Vec_XstripQ[i], Vec_XstripQ[i + 1]);
2906 swap(Vec_VstripQ[i], Vec_VstripQ[i + 1]);
2907 swap(Vec_XstripT[i], Vec_XstripT[i + 1]);
2908 swap(Vec_VstripT[i], Vec_VstripT[i + 1]);
2909 swap(Vec_XstripTf[i], Vec_XstripTf[i + 1]);
2910 swap(Vec_VstripTf[i], Vec_VstripTf[i + 1]);
2911 swap(Vec_XstripID[i], Vec_XstripID[i + 1]);
2912 swap(Vec_VstripID[i], Vec_VstripID[i + 1]);
2913 } else if (Vec_flag[i] == Vec_flag[i + 1] && Vec_layerid[i] == Vec_layerid[i + 1] && (Vec_phi[i] < Vec_phi[i + 1])) {
2914 Swap(i, i + 1);
2915 swap(Vec_xCluSize[i], Vec_xCluSize[i + 1]);
2916 swap(Vec_vCluSize[i], Vec_vCluSize[i + 1]);
2917 swap(Vec_XQ[i], Vec_XQ[i + 1]);
2918 swap(Vec_VQ[i], Vec_VQ[i + 1]);
2919 swap(Vec_XstripQ[i], Vec_XstripQ[i + 1]);
2920 swap(Vec_VstripQ[i], Vec_VstripQ[i + 1]);
2921 swap(Vec_XstripT[i], Vec_XstripT[i + 1]);
2922 swap(Vec_VstripT[i], Vec_VstripT[i + 1]);
2923 swap(Vec_XstripTf[i], Vec_XstripTf[i + 1]);
2924 swap(Vec_VstripTf[i], Vec_VstripTf[i + 1]);
2925 swap(Vec_XstripID[i], Vec_XstripID[i + 1]);
2926 swap(Vec_VstripID[i], Vec_VstripID[i + 1]);
2927 } else {
2928 noswap++;
2929 }
2930 }
2931 }
2932}
2933
2934bool CgemLineFit::Erase_outer() {
2935
2936 vector<double> _x, _v;
2937 vector<int> _iter;
2938 int wrong(-1);
2939 for (int i = 0; i < 3; i++) {
2940 _x.push_back(Vec_phi[i]);
2941 _v.push_back(Vec_v[i]);
2942 _iter.push_back(i);
2943 }
2944 if (_x[0] == _x[1] && _v[0] == _v[2]) wrong = 0;
2945 else if (_x[0] == _x[1] && _v[1] == _v[2]) wrong = 1;
2946 else if (_x[0] == _x[2] && _v[0] == _v[1]) wrong = 0;
2947 else if (_x[0] == _x[2] && _v[2] == _v[1]) wrong = 2;
2948 else if (_x[1] == _x[2] && _v[1] == _v[0]) wrong = 1;
2949 else if (_x[1] == _x[2] && _v[2] == _v[0]) wrong = 2;
2950 else { return false; }
2951 erase(_iter[wrong]);
2952 return true;
2953}
2954
2955void CgemLineFit::erase(int i) {
2956 vector<double>::iterator _iter_v;
2957 _iter_v = Vec_v.begin();
2958 Vec_v.erase(_iter_v + i);
2959
2960 vector<double>::iterator _iter_phi;
2961 _iter_phi = Vec_phi.begin();
2962 Vec_phi.erase(_iter_phi + i);
2963
2964 vector<double>::iterator _iter_Z;
2965 _iter_Z = Vec_Z.begin();
2966 Vec_Z.erase(_iter_Z + i);
2967
2968 vector<double>::iterator _iter_layer;
2969 _iter_layer = Vec_layer.begin();
2970 Vec_layer.erase(_iter_layer + i);
2971
2972 vector<double>::iterator _iter_layerid;
2973 _iter_layerid = Vec_layerid.begin();
2974 Vec_layerid.erase(_iter_layerid + i);
2975
2976 vector<double>::iterator _iter_Virlayerid;
2977 _iter_Virlayerid = Vec_Virlayerid.begin();
2978 Vec_Virlayerid.erase(_iter_Virlayerid + i);
2979
2980 vector<double>::iterator _iter_flag;
2981 _iter_flag = Vec_flag.begin();
2982 Vec_flag.erase(_iter_flag + i);
2983}
2984
2985void CgemLineFit::Discard(int layer) {
2986 for (int i = 0; i < Vec_layer.size(); i++) {
2987 if (Vec_layerid[i] == layer) erase(i);
2988 }
2989}
2990
2991StraightLine *CgemLineFit::IniPar(double phi1, double z1, int i, double phi2, double z2, int j) {
2992 int Gi = int(i / 2);
2993 int Gj = int(j / 2);
2994
2995 HepPoint3D up(R_layer[Gi] * cos(phi1), R_layer[Gi] * sin(phi1), z1);
2996 HepPoint3D down(R_layer[Gj] * cos(phi2), R_layer[Gj] * sin(phi2), z2);
2997
2998 if (Align_Flag) up = Al->point_invTransform(i, up); // added by Guoaq-2020/03/13
2999 if (Align_Flag) down = Al->point_invTransform(j, down);
3000
3001 StraightLine *l = new StraightLine(up, down);
3002 return l;
3003}
3004
3005// I guess this function is designed to remove the mis-combination of the 2D clusters on the layer0
3006// I think the idea is to compare the distance of clusters to the intersections of the given line, and chose the pair of closest
3007// one
3008void CgemLineFit::Filter(int layerid, StraightLine *l1) {
3009 if (!Layer_cross(layerid, l1)) { Discard(layerid); }
3010
3011 vector<int> number;
3012 vector<double> dst;
3013 vector<double> dst_up;
3014 vector<double> dst_down;
3015 vector<double> v_x;
3016 vector<double> v_v;
3017
3018 if (Run10_Flag) {
3019 for (int i = 0; i < Vec_layer.size(); i++) {
3020 if (Vec_layerid[i] != layerid) continue;
3021 number.push_back(i);
3022 // obtain the distance between the up/down side of crosssections to the given position (phi, V)
3023 dst_up.push_back(Min_dis_up(layerid, Vec_phi[i], Vec_v[i], l1));
3024 dst_down.push_back(Min_dis_down(layerid, Vec_phi[i], Vec_v[i], l1));
3025 }
3026
3027 for (int j = 0; j < number.size(); j++) { int i = number[j]; }
3028
3029 // sort the clusters according to the up side distance
3030 int Loop;
3031 Loop = 0;
3032 while (Loop != number.size() - 1) {
3033 Loop = 0;
3034 for (int j = 0; j < number.size() - 1; j++) {
3035 if (dst_up[j] > dst_up[j + 1]) {
3036 swap(dst_up[j], dst_up[j + 1]);
3037 swap(dst_down[j], dst_down[j + 1]);
3038 swap(number[j], number[j + 1]);
3039 } else Loop++;
3040 }
3041 }
3042
3043 // set the index of the nearest cluster
3044 for (int j = 0; j < number.size(); j++) { int i = number[j]; }
3045
3046 // remove the first member
3047 number.erase(number.begin());
3048 dst_up.erase(dst_up.begin());
3049 dst_down.erase(dst_down.begin());
3050
3051 // what's the difference to previous operation?
3052 for (int j = 0; j < number.size(); j++) { int i = number[j]; }
3053
3054 // sort the clusters according to the down side distance
3055 Loop = 0;
3056 while (Loop != number.size() - 1) {
3057 Loop = 0;
3058 for (int j = 0; j < number.size() - 1; j++) {
3059 if (dst_down[j] > dst_down[j + 1]) {
3060 swap(dst_down[j], dst_down[j + 1]);
3061 swap(number[j], number[j + 1]);
3062 } else Loop++;
3063 }
3064 }
3065
3066 for (int j = 0; j < number.size(); j++) { int i = number[j]; }
3067
3068 number.erase(number.begin());
3069
3070 } else // if run10_flag is false
3071 {
3072 for (int i = 0; i < Vec_layer.size(); i++) {
3073 if (Vec_layerid[i] != layerid) continue;
3074 number.push_back(i);
3075
3076 dst.push_back(Min_dis(layerid, Vec_phi[i], Vec_v[i], l1));
3077 }
3078
3079 int Loop(0);
3080
3081 while (Loop != number.size() - 1) {
3082 Loop = 0;
3083 for (int j = 0; j < number.size() - 1; j++) {
3084 if (dst[j] > dst[j + 1]) {
3085 swap(dst[j], dst[j + 1]);
3086 swap(number[j], number[j + 1]);
3087 } else Loop++;
3088 }
3089 }
3090
3091 number.erase(number.begin());
3092
3093 } // end of if run10_flag
3094
3095 for (int j = 0; j < number.size(); j++) { int i = number[j]; }
3096
3097 sort(number.begin(), number.end(), less<int>());
3098 for (int i = number.size() - 1; i != -1; i--) //{ erase(number[i]); }//removed by Farinelli Riccardo in 2023.12.29
3099
3100 return;
3101}
3102
3103bool CgemLineFit::Get_MCInf() {
3104 vector<double> Vec_truth;
3105 Vec_truth.clear();
3106 Vec_truth = MC_truth();
3107
3108 vector<double> _Vec_truth = Trans(Vec_truth[0], Vec_truth[1], Vec_truth[2], Vec_truth[3]);
3109
3110 MC_par0 = _Vec_truth[0];
3111 MC_par1 = _Vec_truth[1];
3112 MC_par2 = _Vec_truth[2];
3113 MC_par3 = _Vec_truth[3];
3114 MC_px = Vec_truth[4];
3115 MC_py = Vec_truth[5];
3116 MC_pz = Vec_truth[6];
3117
3118 vector<double> V_MC_clusters = Get_Clusters(Vec_truth);
3119 x_0_up_mc = V_MC_clusters[0];
3120 v_0_up_mc = V_MC_clusters[1];
3121 x_0_down_mc = V_MC_clusters[2];
3122 v_0_down_mc = V_MC_clusters[3];
3123
3124 x_1_up_mc = V_MC_clusters[4];
3125 v_1_up_mc = V_MC_clusters[5];
3126 x_1_down_mc = V_MC_clusters[6];
3127 v_1_down_mc = V_MC_clusters[7];
3128
3129 x_2_up_mc = V_MC_clusters[8];
3130 v_2_up_mc = V_MC_clusters[9];
3131 x_2_down_mc = V_MC_clusters[10];
3132 v_2_down_mc = V_MC_clusters[11];
3133 return true;
3134}
3135
3136vector<double> CgemLineFit::MC_truth() {
3137 int nTruth[3] = {0, 0, 0};
3138 int trkID_Layer1[100], trkID_Layer2[100], trkID_Layer3[100];
3139 int pdg_Layer1[100], pdg_Layer2[100], pdg_Layer3[100];
3140 double x_pre_Layer1[100], x_pre_Layer2[100], x_pre_Layer3[100];
3141 double y_pre_Layer1[100], y_pre_Layer2[100], y_pre_Layer3[100];
3142 double z_pre_Layer1[100], z_pre_Layer2[100], z_pre_Layer3[100];
3143 double x_post_Layer1[100], x_post_Layer2[100], x_post_Layer3[100];
3144 double y_post_Layer1[100], y_post_Layer2[100], y_post_Layer3[100];
3145 double z_post_Layer1[100], z_post_Layer2[100], z_post_Layer3[100];
3146 double phi_Layer1[100], phi_Layer2[100], phi_Layer3[100];
3147 double z_Layer1[100], z_Layer2[100], z_Layer3[100];
3148 double V_Layer1[100], V_Layer2[100], V_Layer3[100];
3149 double px_pre_Layer1[100], px_pre_Layer2[100], px_pre_Layer3[100];
3150 double py_pre_Layer1[100], py_pre_Layer2[100], py_pre_Layer3[100];
3151 double pz_pre_Layer1[100], pz_pre_Layer2[100], pz_pre_Layer3[100];
3152 double en_Layer1[100], en_Layer2[100], en_Layer3[100];
3153 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
3154 /* if (!cgemMcHitCol)
3155 {
3156 std::cout << "Could not retrieve cgemMcHitCol" << std::endl;
3157 return;
3158 } */
3159 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
3160 for (; iter_truth != cgemMcHitCol->end(); iter_truth++) {
3161 int iLayer = (*iter_truth)->GetLayerID();
3162 int trkID = (*iter_truth)->GetTrackID();
3163 int pdg = (*iter_truth)->GetPDGCode();
3164 double x_pre = (*iter_truth)->GetPositionXOfPrePoint(); // in mm
3165 double y_pre = (*iter_truth)->GetPositionYOfPrePoint();
3166 double z_pre = (*iter_truth)->GetPositionZOfPrePoint();
3167 double x_post = (*iter_truth)->GetPositionXOfPostPoint();
3168 double y_post = (*iter_truth)->GetPositionYOfPostPoint();
3169 double z_post = (*iter_truth)->GetPositionZOfPostPoint();
3170 double x_middle = 0.5 * (x_pre + x_post);
3171 double y_middle = 0.5 * (y_pre + y_post);
3172 double z_middle = 0.5 * (z_pre + z_post);
3173 double r_middle = sqrt(x_middle * x_middle + y_middle * y_middle); // cout<<"r_middle = "<<r_middle<<endl;
3174 double phi_middle = atan2(y_middle, x_middle);
3175 double px_pre = (*iter_truth)->GetMomentumXOfPrePoint();
3176 double py_pre = (*iter_truth)->GetMomentumYOfPrePoint();
3177 double pz_pre = (*iter_truth)->GetMomentumZOfPrePoint();
3178 double en = (*iter_truth)->GetTotalEnergyDeposit();
3179 CgemGeoReadoutPlane *readoutPlane = NULL;
3180 // cout<<"phi_middle, z_middle = "<<phi_middle<<", "<<z_middle<<endl;
3181 // cout<<"iLayer = "<<iLayer<<endl;
3182 for (int j = 0; j < myNSheets[iLayer]; j++) {
3183 // cout<<"j = "<<j <<endl;
3184 readoutPlane = m_cgemGeomSvc->getReadoutPlane(iLayer, j);
3185 // cout<<"OnThePlane: "<<readoutPlane->OnThePlane(phi_middle,z_middle)<<endl;
3186 if (readoutPlane->OnThePlane(phi_middle, z_middle)) break;
3187 readoutPlane = NULL;
3188 }
3189 double V_mid = 9999;
3190 if (readoutPlane == NULL) {
3191 cout << "CgemClusterCreate::fillMCTruth() readoutPlane not found with phi_middle = " << phi_middle << ", layer = " << iLayer
3192 << endl;
3193 } else {
3194 // readoutPlane->print();
3195 // cout<<"phi_middle, z_middle = "<<phi_middle<<", "<<z_middle<<endl;
3196 V_mid = readoutPlane->getVFromPhiZ(phi_middle, z_middle);
3197 // cout<<"V_mid = "<<V_mid<<endl;
3198 // double V_mid = readoutPlane->GetVFromLocalXZ(phi_middle,z_middle);
3199 }
3200 switch (iLayer) {
3201 case 0:
3202 if (nTruth[0] >= 100) break;
3203 trkID_Layer1[nTruth[0]] = trkID;
3204 pdg_Layer1[nTruth[0]] = pdg;
3205 x_pre_Layer1[nTruth[0]] = x_pre;
3206 y_pre_Layer1[nTruth[0]] = y_pre;
3207 z_pre_Layer1[nTruth[0]] = z_pre;
3208 x_post_Layer1[nTruth[0]] = x_post;
3209 y_post_Layer1[nTruth[0]] = y_post;
3210 z_post_Layer1[nTruth[0]] = z_post;
3211 phi_Layer1[nTruth[0]] = atan2(y_post + y_pre, x_post + x_pre);
3212 z_Layer1[nTruth[0]] = z_middle;
3213 V_Layer1[nTruth[0]] = V_mid;
3214 if (atan2(y_post + y_pre, x_post + x_pre) > 0) {
3215 x_0_up_mc2 = atan2(y_post + y_pre, x_post + x_pre);
3216 z_0_up_mc2 = z_middle;
3217 v_0_up_mc2 = V_mid;
3218 } else {
3219 x_0_down_mc2 = atan2(y_post + y_pre, x_post + x_pre);
3220 z_0_down_mc2 = z_middle;
3221 v_0_down_mc2 = V_mid;
3222 }
3223 px_pre_Layer1[nTruth[0]] = px_pre;
3224 py_pre_Layer1[nTruth[0]] = py_pre;
3225 pz_pre_Layer1[nTruth[0]] = pz_pre;
3226 en_Layer1[nTruth[0]] = en;
3227 break;
3228 case 1:
3229 if (nTruth[1] >= 100) break;
3230 trkID_Layer2[nTruth[1]] = trkID;
3231 pdg_Layer2[nTruth[1]] = pdg;
3232 x_pre_Layer2[nTruth[1]] = x_pre;
3233 y_pre_Layer2[nTruth[1]] = y_pre;
3234 z_pre_Layer2[nTruth[1]] = z_pre;
3235 x_post_Layer2[nTruth[1]] = x_post;
3236 y_post_Layer2[nTruth[1]] = y_post;
3237 z_post_Layer2[nTruth[1]] = z_post;
3238 phi_Layer2[nTruth[1]] = atan2(y_post + y_pre, x_post + x_pre);
3239 z_Layer2[nTruth[1]] = z_middle;
3240 V_Layer2[nTruth[1]] = V_mid;
3241 if (atan2(y_post + y_pre, x_post + x_pre) > 0) {
3242 x_1_up_mc2 = atan2(y_post + y_pre, x_post + x_pre);
3243 z_1_up_mc2 = z_middle;
3244 v_1_up_mc2 = V_mid;
3245 } else {
3246 x_1_down_mc2 = atan2(y_post + y_pre, x_post + x_pre);
3247 z_1_down_mc2 = z_middle;
3248 v_1_down_mc2 = V_mid;
3249 }
3250 px_pre_Layer2[nTruth[1]] = px_pre;
3251 py_pre_Layer2[nTruth[1]] = py_pre;
3252 pz_pre_Layer2[nTruth[1]] = pz_pre;
3253 en_Layer2[nTruth[1]] = en;
3254 break;
3255 case 2:
3256 if (nTruth[2] >= 100) break;
3257 trkID_Layer3[nTruth[2]] = trkID;
3258 pdg_Layer3[nTruth[2]] = pdg;
3259 x_pre_Layer3[nTruth[2]] = x_pre;
3260 y_pre_Layer3[nTruth[2]] = y_pre;
3261 z_pre_Layer3[nTruth[2]] = z_pre;
3262 x_post_Layer3[nTruth[2]] = x_post;
3263 y_post_Layer3[nTruth[2]] = y_post;
3264 z_post_Layer3[nTruth[2]] = z_post;
3265 phi_Layer3[nTruth[2]] = atan2(y_post + y_pre, x_post + x_pre);
3266 z_Layer3[nTruth[2]] = z_middle;
3267 V_Layer3[nTruth[2]] = V_mid;
3268 if (atan2(y_post + y_pre, x_post + x_pre) > 0) {
3269 x_2_up_mc2 = atan2(y_post + y_pre, x_post + x_pre);
3270 v_2_up_mc2 = V_mid;
3271 z_2_up_mc2 = z_middle;
3272 } else {
3273 x_2_down_mc2 = atan2(y_post + y_pre, x_post + x_pre);
3274 v_2_down_mc2 = V_mid;
3275 z_2_down_mc2 = z_middle;
3276 }
3277 px_pre_Layer3[nTruth[2]] = px_pre;
3278 py_pre_Layer3[nTruth[2]] = py_pre;
3279 pz_pre_Layer3[nTruth[2]] = pz_pre;
3280 en_Layer3[nTruth[2]] = en;
3281 break;
3282 default: cout << "wrong layer ID for CGEM: " << iLayer << endl;
3283 }
3284 nTruth[iLayer]++;
3285 } // loop hit truth
3286 for (int i = 0; i < 3; i++)
3287 if (nTruth[i] > 100) nTruth[i] = 100;
3288 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
3289
3290 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
3291 HepLorentzVector p4_mu(0, 0, 0, 0);
3292 HepLorentzVector p4_pos(0, 0, 0, 0);
3293 for (; iter_mc != mcParticleCol->end(); iter_mc++) {
3294 if (!(*iter_mc)->decayFromGenerator()) continue;
3295 HepLorentzVector p4_mc(0, 0, 0, 0);
3296
3297 if (13 == abs((*iter_mc)->particleProperty())) {
3298
3299 p4_mu = (*iter_mc)->initialFourMomentum();
3300 p4_pos = (*iter_mc)->initialPosition();
3301 }
3302 }
3303 vector<double> _mc;
3304 Hep3Vector _BP(p4_pos.x() * 10, p4_pos.y() * 10, p4_pos.z() * 10);
3305 pos_x = _BP[0];
3306 pos_y = _BP[1];
3307 pos_z = _BP[2];
3308 _mc = Get4Par(p4_mu, _BP);
3309
3310 _mc.push_back(p4_mu.px());
3311 _mc.push_back(p4_mu.py());
3312 _mc.push_back(p4_mu.pz());
3313 if (HitPos(p4_mu, _BP)) _mc.push_back(1);
3314 else _mc.push_back(0);
3315
3316 return _mc;
3317}
3318
3319vector<double> CgemLineFit::Get4Par(HepLorentzVector p4, Hep3Vector bp) {
3320 vector<double> _V_par;
3321 _V_par.clear();
3322 double xb = p4.px();
3323 double yb = p4.py();
3324 double zb = p4.pz();
3325 double xa = bp[0];
3326 double ya = bp[1];
3327 double za = bp[2];
3328 double k = -1 * (xa * xb + ya * yb) / (xb * xb + yb * yb);
3329 double xc = k * xb + xa;
3330 double yc = k * yb + ya;
3331 double zc = k * zb + za;
3332 HepPoint3D a(xa, ya, za);
3333 HepPoint3D c(xc, yc, zc);
3334 StraightLine l(a, c);
3335 _V_par.push_back(l.dr());
3336 _V_par.push_back(l.phi0());
3337 _V_par.push_back(l.dz());
3338 _V_par.push_back(l.tanl());
3339 return _V_par;
3340}
3341
3342bool CgemLineFit::HitPos(HepLorentzVector p4, Hep3Vector bp) {
3343 double L(500), W(160), H(600);
3344
3345 vector<double> _V_par;
3346 _V_par.clear();
3347 double xb = p4.px();
3348 double yb = p4.py();
3349 double zb = p4.pz();
3350 double xa = bp[0];
3351 double ya = bp[1];
3352 double za = bp[2];
3353 double k = H / yb * -1;
3354 double xc = k * xb + xa;
3355 double yc = k * yb + ya;
3356 double zc = k * zb + za;
3357 hit_x = xc;
3358 hit_y = yc;
3359 hit_z = zc;
3360
3361 if (fabs(xc) > W / 2) return false;
3362 if (fabs(zc) > L / 2) return false;
3363
3364 else return true;
3365}
3366
3367void CgemLineFit::Get_OtherIniPar(int clst_0, int clst_1, int clst_2) {
3368
3369 if (clst_2 == 0) { // only 2 layers situation
3370
3371 if (clst_0 >= 2 && clst_1 >= 2) {
3372 StraightLine *l_ini_i = IniPar(Vec_phi[2], Vec_Z[2], 1, Vec_phi[3], Vec_Z[3],
3373 0); // why use index of 4 and 5 because we only have 2 layer. I guess hongpeng want to get
3374 // the initial parameters by the innermost layer
3375
3376 vector<double> inii_par = Trans(l_ini_i->dr(), l_ini_i->phi0(), l_ini_i->dz(), l_ini_i->tanl());
3377 inii_0 = inii_par[0];
3378 inii_1 = inii_par[1];
3379 inii_2 = inii_par[2];
3380 inii_3 = inii_par[3];
3381 delete l_ini_i;
3382 }
3383 if (clst_1 >= 2) {
3384 StraightLine *l_ini_m = IniPar(Vec_phi[0], Vec_Z[0], 3, Vec_phi[1], Vec_Z[1],
3385 2); // I guess Hongpeng want to get the initial parameters by the middle layer
3386 vector<double> inim_par = Trans(l_ini_m->dr(), l_ini_m->phi0(), l_ini_m->dz(), l_ini_m->tanl());
3387
3388 inim_0 = inim_par[0];
3389 inim_1 = inim_par[1];
3390 inim_2 = inim_par[2];
3391 inim_3 = inim_par[3];
3392 delete l_ini_m;
3393 }
3394
3395 } else {
3396 if (clst_0 >= 2 && clst_1 >= 2) {
3397 StraightLine *l_ini_i = IniPar(Vec_phi[4], Vec_Z[4], 1, Vec_phi[5], Vec_Z[5],
3398 0); // why use index of 4 and 5 because we only have 2 layer. I guess hongpeng want to get
3399 // the initial parameters by the innermost layer
3400
3401 vector<double> inii_par = Trans(l_ini_i->dr(), l_ini_i->phi0(), l_ini_i->dz(), l_ini_i->tanl());
3402 inii_0 = inii_par[0];
3403 inii_1 = inii_par[1];
3404 inii_2 = inii_par[2];
3405 inii_3 = inii_par[3];
3406 delete l_ini_i;
3407 }
3408 if (clst_1 >= 2) {
3409 StraightLine *l_ini_m = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3],
3410 2); // I guess Hongpeng want to get the initial parameters by the middle layer
3411 vector<double> inim_par = Trans(l_ini_m->dr(), l_ini_m->phi0(), l_ini_m->dz(), l_ini_m->tanl());
3412
3413 inim_0 = inim_par[0];
3414 inim_1 = inim_par[1];
3415 inim_2 = inim_par[2];
3416 inim_3 = inim_par[3];
3417 delete l_ini_m;
3418 }
3419 }
3420}
3421
3422void CgemLineFit::Rec_Clusters() {
3423
3424 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
3425
3426 if (Vec_Virlayerid[i] == 5) {
3427 x_2_up = Vec_phi[i];
3428 v_2_up = Vec_v[i];
3429 z_2_up = Vec_Z[i];
3430 }
3431 if (Vec_Virlayerid[i] == 4) {
3432 x_2_down = Vec_phi[i];
3433 v_2_down = Vec_v[i];
3434 z_2_down = Vec_Z[i];
3435 }
3436 if (Vec_Virlayerid[i] == 3) {
3437 x_1_up = Vec_phi[i];
3438 v_1_up = Vec_v[i];
3439 z_1_up = Vec_Z[i];
3440 }
3441
3442 if (Vec_Virlayerid[i] == 2) {
3443 x_1_down = Vec_phi[i];
3444 v_1_down = Vec_v[i];
3445 z_1_down = Vec_Z[i];
3446 }
3447
3448 if (Vec_Virlayerid[i] == 1) {
3449 x_0_up = Vec_phi[i];
3450 v_0_up = Vec_v[i];
3451 z_0_up = Vec_Z[i];
3452 }
3453
3454 if (Vec_Virlayerid[i] == 0) {
3455 x_0_down = Vec_phi[i];
3456 v_0_down = Vec_v[i];
3457 z_0_down = Vec_Z[i];
3458 }
3459 }
3460}
3461
3462void CgemLineFit::Rec_Clusters_mTPC() {
3463
3464 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
3465 if (Vec_Virlayerid[i] == 5) {
3466 x_2_up_m = Vec_m_phi[i];
3467 v_2_up_m = Vec_m_v[i];
3468 z_2_up_m = Vec_m_Z[i];
3469 }
3470 if (Vec_Virlayerid[i] == 4) {
3471 x_2_down_m = Vec_m_phi[i];
3472 v_2_down_m = Vec_m_v[i];
3473 z_2_down_m = Vec_m_Z[i];
3474 }
3475 if (Vec_Virlayerid[i] == 3) {
3476 x_1_up_m = Vec_m_phi[i];
3477 v_1_up_m = Vec_m_v[i];
3478 z_1_up_m = Vec_m_Z[i];
3479 }
3480
3481 if (Vec_Virlayerid[i] == 2) {
3482 x_1_down_m = Vec_m_phi[i];
3483 v_1_down_m = Vec_m_v[i];
3484 z_1_down_m = Vec_m_Z[i];
3485 }
3486
3487 if (Vec_Virlayerid[i] == 1) {
3488 x_0_up_m = Vec_m_phi[i];
3489 v_0_up_m = Vec_m_v[i];
3490 z_0_up_m = Vec_m_Z[i];
3491 }
3492 if (Vec_Virlayerid[i] == 0) {
3493 x_0_down_m = Vec_m_phi[i];
3494 v_0_down_m = Vec_m_v[i];
3495 z_0_down_m = Vec_m_Z[i];
3496 }
3497 }
3498}
3499void CgemLineFit::Rec_ClusterSize() {
3500 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
3501 if (Vec_Virlayerid[i] == 5) {
3502 x_2_up_size = Vec_xCluSize[i];
3503 v_2_up_size = Vec_vCluSize[i];
3504 }
3505 if (Vec_Virlayerid[i] == 4) {
3506 x_2_down_size = Vec_xCluSize[i];
3507 v_2_down_size = Vec_vCluSize[i];
3508 }
3509 if (Vec_Virlayerid[i] == 3) {
3510 x_1_up_size = Vec_xCluSize[i];
3511 v_1_up_size = Vec_vCluSize[i];
3512 }
3513 if (Vec_Virlayerid[i] == 2) {
3514 x_1_down_size = Vec_xCluSize[i];
3515 v_1_down_size = Vec_vCluSize[i];
3516 }
3517 if (Vec_Virlayerid[i] == 1) {
3518 x_0_up_size = Vec_xCluSize[i];
3519 v_0_up_size = Vec_vCluSize[i];
3520 }
3521 if (Vec_Virlayerid[i] == 0) {
3522 x_0_down_size = Vec_xCluSize[i];
3523 v_0_down_size = Vec_vCluSize[i];
3524 }
3525 }
3526}
3527
3528void CgemLineFit::Rec_ClusterQ() {
3529 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
3530 if (Vec_Virlayerid[i] == 5) {
3531 x_2_up_Q = Vec_XQ[i];
3532 v_2_up_Q = Vec_VQ[i];
3533 }
3534 if (Vec_Virlayerid[i] == 4) {
3535 x_2_down_Q = Vec_XQ[i];
3536 v_2_down_Q = Vec_VQ[i];
3537 }
3538 if (Vec_Virlayerid[i] == 3) {
3539 x_1_up_Q = Vec_XQ[i];
3540 v_1_up_Q = Vec_VQ[i];
3541 }
3542 if (Vec_Virlayerid[i] == 2) {
3543 x_1_down_Q = Vec_XQ[i];
3544 v_1_down_Q = Vec_VQ[i];
3545 }
3546 if (Vec_Virlayerid[i] == 1) {
3547 x_0_up_Q = Vec_XQ[i];
3548 v_0_up_Q = Vec_VQ[i];
3549 }
3550 if (Vec_Virlayerid[i] == 0) {
3551 x_0_down_Q = Vec_XQ[i];
3552 v_0_down_Q = Vec_VQ[i];
3553 }
3554 }
3555}
3556
3557void CgemLineFit::fireStripQ() {
3558 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
3559 if (Vec_Virlayerid[i] == 5) {
3560 for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_2_up_stripQ[j] = Vec_XstripQ[i][j]; }
3561 for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_2_up_stripQ[j] = Vec_VstripQ[i][j]; }
3562 }
3563 if (Vec_Virlayerid[i] == 4) {
3564 for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_2_down_stripQ[j] = Vec_XstripQ[i][j]; }
3565 for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_2_down_stripQ[j] = Vec_VstripQ[i][j]; }
3566 }
3567 if (Vec_Virlayerid[i] == 3) {
3568 for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_1_up_stripQ[j] = Vec_XstripQ[i][j]; }
3569 for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_1_up_stripQ[j] = Vec_VstripQ[i][j]; }
3570 }
3571 if (Vec_Virlayerid[i] == 2) {
3572 for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_1_down_stripQ[j] = Vec_XstripQ[i][j]; }
3573 for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_1_down_stripQ[j] = Vec_VstripQ[i][j]; }
3574 }
3575 if (Vec_Virlayerid[i] == 1) {
3576 for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_0_up_stripQ[j] = Vec_XstripQ[i][j]; }
3577 for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_0_up_stripQ[j] = Vec_VstripQ[i][j]; }
3578 }
3579 if (Vec_Virlayerid[i] == 0) {
3580 for (int j = 0; j < Vec_XstripQ[i].size(); j++) { x_0_down_stripQ[j] = Vec_XstripQ[i][j]; }
3581 for (int j = 0; j < Vec_VstripQ[i].size(); j++) { v_0_down_stripQ[j] = Vec_VstripQ[i][j]; }
3582 }
3583 }
3584}
3585
3586void CgemLineFit::fireStripT() {
3587 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
3588 if (Vec_Virlayerid[i] == 5) {
3589 for (int j = 0; j < Vec_XstripT[i].size(); j++) {
3590 x_2_up_stripT[j] = Vec_XstripT[i][j];
3591 x_2_up_stripTf[j] = Vec_XstripTf[i][j];
3592 }
3593 for (int j = 0; j < Vec_VstripT[i].size(); j++) {
3594 v_2_up_stripT[j] = Vec_VstripT[i][j];
3595 v_2_up_stripTf[j] = Vec_VstripTf[i][j];
3596 }
3597 }
3598 if (Vec_Virlayerid[i] == 4) {
3599 for (int j = 0; j < Vec_XstripT[i].size(); j++) {
3600 x_2_down_stripT[j] = Vec_XstripT[i][j];
3601 x_2_down_stripTf[j] = Vec_XstripTf[i][j];
3602 }
3603 for (int j = 0; j < Vec_VstripT[i].size(); j++) {
3604 v_2_down_stripT[j] = Vec_VstripT[i][j];
3605 v_2_down_stripTf[j] = Vec_VstripTf[i][j];
3606 }
3607 }
3608 if (Vec_Virlayerid[i] == 3) {
3609 for (int j = 0; j < Vec_XstripT[i].size(); j++) {
3610 x_1_up_stripT[j] = Vec_XstripT[i][j];
3611 x_1_up_stripTf[j] = Vec_XstripTf[i][j];
3612 }
3613 for (int j = 0; j < Vec_VstripT[i].size(); j++) {
3614 v_1_up_stripT[j] = Vec_VstripT[i][j];
3615 v_1_up_stripTf[j] = Vec_VstripTf[i][j];
3616 }
3617 }
3618 if (Vec_Virlayerid[i] == 2) {
3619 for (int j = 0; j < Vec_XstripT[i].size(); j++) {
3620 x_1_down_stripT[j] = Vec_XstripT[i][j];
3621 x_1_down_stripTf[j] = Vec_XstripTf[i][j];
3622 }
3623 for (int j = 0; j < Vec_VstripT[i].size(); j++) {
3624 v_1_down_stripT[j] = Vec_VstripT[i][j];
3625 v_1_down_stripTf[j] = Vec_VstripTf[i][j];
3626 }
3627 }
3628 if (Vec_Virlayerid[i] == 1) {
3629 for (int j = 0; j < Vec_XstripT[i].size(); j++) {
3630 x_0_up_stripT[j] = Vec_XstripT[i][j];
3631 x_0_up_stripTf[j] = Vec_XstripTf[i][j];
3632 }
3633 for (int j = 0; j < Vec_VstripT[i].size(); j++) {
3634 v_0_up_stripT[j] = Vec_VstripT[i][j];
3635 v_0_up_stripTf[j] = Vec_VstripTf[i][j];
3636 }
3637 }
3638 if (Vec_Virlayerid[i] == 0) {
3639 for (int j = 0; j < Vec_XstripT[i].size(); j++) {
3640 x_0_down_stripT[j] = Vec_XstripT[i][j];
3641 x_0_down_stripTf[j] = Vec_XstripTf[i][j];
3642 }
3643 for (int j = 0; j < Vec_VstripT[i].size(); j++) {
3644 v_0_down_stripT[j] = Vec_VstripT[i][j];
3645 v_0_down_stripTf[j] = Vec_VstripTf[i][j];
3646 }
3647 }
3648 }
3649}
3650
3651void CgemLineFit::fireStripID() {
3652 for (int i = 0; i < Vec_Virlayerid.size(); i++) {
3653 if (Vec_Virlayerid[i] == 5) {
3654 for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_2_up_stripID[j] = Vec_XstripID[i][j]; }
3655 for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_2_up_stripID[j] = Vec_VstripID[i][j]; }
3656 }
3657 if (Vec_Virlayerid[i] == 4) {
3658 for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_2_down_stripID[j] = Vec_XstripID[i][j]; }
3659 for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_2_down_stripID[j] = Vec_VstripID[i][j]; }
3660 }
3661 if (Vec_Virlayerid[i] == 3) {
3662 for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_1_up_stripID[j] = Vec_XstripID[i][j]; }
3663 for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_1_up_stripID[j] = Vec_VstripID[i][j]; }
3664 }
3665 if (Vec_Virlayerid[i] == 2) {
3666 for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_1_down_stripID[j] = Vec_XstripID[i][j]; }
3667 for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_1_down_stripID[j] = Vec_VstripID[i][j]; }
3668 }
3669 if (Vec_Virlayerid[i] == 1) {
3670 for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_0_up_stripID[j] = Vec_XstripID[i][j]; }
3671 for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_0_up_stripID[j] = Vec_VstripID[i][j]; }
3672 }
3673 if (Vec_Virlayerid[i] == 0) {
3674 for (int j = 0; j < Vec_XstripID[i].size(); j++) { x_0_down_stripID[j] = Vec_XstripID[i][j]; }
3675 for (int j = 0; j < Vec_VstripID[i].size(); j++) { v_0_down_stripID[j] = Vec_VstripID[i][j]; }
3676 }
3677 }
3678}
3679
3680void CgemLineFit::Fit_Clusters(double par[]) {
3681 vector<double> Vec_truth;
3682 Vec_truth.push_back(par[0]);
3683 Vec_truth.push_back(par[1]);
3684 Vec_truth.push_back(par[2]);
3685 Vec_truth.push_back(par[3]);
3686
3687 vector<double> clusters = Get_Clusters(Vec_truth);
3688
3689 x_0_up_f = clusters[0];
3690 v_0_up_f = clusters[1];
3691 x_0_down_f = clusters[2];
3692 v_0_down_f = clusters[3];
3693
3694 x_1_up_f = clusters[4];
3695 v_1_up_f = clusters[5];
3696 x_1_down_f = clusters[6];
3697 v_1_down_f = clusters[7];
3698
3699 x_2_up_f = clusters[8];
3700 v_2_up_f = clusters[9];
3701 x_2_down_f = clusters[10];
3702 v_2_down_f = clusters[11];
3703}
3704
3705void CgemLineFit::GetIntersection(StraightLine* line) {
3706
3707 if (debug) cout << "Start Getintersection" << endl;
3708 Double_t fmin, edm, errdef;
3709 Int_t ierflg, npari, nparx, istat;
3710 double trkpar[4] = {-999, -999, -999, -999};
3711 double trkparErr[4] = {-999, -999, -999, -999};
3712
3713 double phiV_up[2], phiV_down[2];
3714 HepPoint3D Up, Down;
3715 double phiV_up_Tmp[2], phiV_down_Tmp[2];
3716 HepPoint3D Up_Tmp, Down_Tmp;
3717
3718 double Ang[3];
3719
3720 //StraightLine *l1a = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
3721 TMinuit *fit = Fit(line, fa - f21, 0);
3722 fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
3723 inter_1_flag = istat;
3724 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
3725 StraightLine *l1b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
3726 Store_Trk(fit, 1000);
3727 delete fit;
3728
3729 if (Align_FLAG) {
3730 Mp->getPointAligned(4, l1b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3731 Mp->getPointAligned(5, l1b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3732 } else Mp->getPointIdealGeom(2 * 2, l1b, Up, Down, phiV_up, phiV_down);
3733
3734 to_0_2pi(phiV_up[0]);
3735 to_0_2pi(phiV_down[0]);
3736 inter_1_x = phiV_up[0];
3737 inter_1_V = phiV_up[1];
3738 inter_1_z = Up[2];
3739
3740 /*inter_4_x=phiV_down[0];
3741 inter_4_V=phiV_down[1];
3742 inter_4_z=Down[2];*/
3743
3744 HepPoint3D l2_up(cos(inter_1_x) * R_layer[2], sin(inter_1_x) * R_layer[2], inter_1_z);
3745 // HepPoint3D l1_down(cos(Vec_phi[1])*R_layer[1],sin(Vec_phi[1])*R_layer[1],Vec_Z[1]);
3746
3747 IncidentAngle(l1b, l2_up, pl_21->getStereoAngle(), Ang);
3748 ang_1_x = Ang[0];
3749 ang_1_v = Ang[1];
3750 ang_1 = Ang[2];
3751
3752 /* IncidentAngle(l1b,l1_down,pl_11->getStereoAngle(),Ang);
3753 ang_4=Ang[0];
3754 ang_4_x=Ang[1];
3755 ang_4_v=Ang[2];*/
3756
3757 //delete l1a;
3758 delete l1b;
3759
3760 //StraightLine *l2a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
3761 fit = Fit(line, fa - f11, 0);
3762 fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
3763 inter_2_flag = istat;
3764 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
3765 Store_Trk(fit, 2000);
3766 delete fit;
3767 StraightLine *l2b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
3768 if (Align_FLAG) {
3769 Mp->getPointAligned(2, l2b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3770 Mp->getPointAligned(3, l2b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3771 } else Mp->getPointIdealGeom(1 * 2, l2b, Up, Down, phiV_up, phiV_down);
3772 to_0_2pi(phiV_up[0]);
3773 inter_2_x = phiV_up[0];
3774 inter_2_V = phiV_up[1];
3775 inter_2_z = Up[2];
3776 HepPoint3D l1_up(cos(inter_2_x) * R_layer[1], sin(inter_2_x) * R_layer[1], inter_2_z);
3777 IncidentAngle(l2b, l1_up, pl_11->getStereoAngle(), Ang);
3778 ang_2_x = Ang[0];
3779 ang_2_v = Ang[1];
3780 ang_2 = Ang[2];
3781 //delete l2a;
3782 delete l2b;
3783
3784 // if(Run10_Flag)
3785 // {
3786 //StraightLine *l3a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
3787 fit = Fit(line, fa - f01, 0);
3788 fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
3789 inter_3_flag = istat;
3790 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
3791 Store_Trk(fit, 3000);
3792 delete fit;
3793 StraightLine *l3b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
3794
3795 if (Align_FLAG) {
3796 Mp->getPointAligned(0, l3b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3797 Mp->getPointAligned(1, l3b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3798 } else Mp->getPointIdealGeom(0 * 2, l3b, Up, Down, phiV_up, phiV_down);
3799 to_0_2pi(phiV_up[0]);
3800 inter_3_x = phiV_up[0];
3801 inter_3_V = phiV_up[1];
3802 inter_3_z = Up[2];
3803
3804 HepPoint3D l0_up(cos(inter_3_x) * R_layer[0], sin(inter_3_x) * R_layer[0], inter_3_z);
3805 IncidentAngle(l3b, l0_up, pl_00->getStereoAngle(), Ang);
3806 ang_3_x = Ang[0];
3807 ang_3_v = Ang[1];
3808 ang_3 = Ang[2];
3809 //delete l3a;
3810 delete l3b;
3811 //}
3812
3813 //StraightLine *l4a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
3814 fit = Fit(line, fa - f00, 0);
3815 fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
3816 inter_4_flag = istat;
3817 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
3818 Store_Trk(fit, 4000);
3819 delete fit;
3820 StraightLine *l4b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
3821 if (Align_FLAG) {
3822 Mp->getPointAligned(0, l4b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3823 Mp->getPointAligned(1, l4b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3824 } else Mp->getPointIdealGeom(0 * 2, l4b, Up, Down, phiV_up, phiV_down);
3825
3826 // to_0_2pi(phiV_up[0]);
3827 to_0_2pi(phiV_down[0]);
3828 inter_4_x = phiV_down[0];
3829 inter_4_V = phiV_down[1];
3830 inter_4_z = Down[2];
3831
3832 HepPoint3D l0_down(cos(inter_4_x) * R_layer[0], sin(inter_4_x) * R_layer[0], inter_4_z);
3833 IncidentAngle(l4b, l0_down, pl_00->getStereoAngle(), Ang);
3834 ang_4_x = Ang[0];
3835 ang_4_v = Ang[1];
3836 ang_4 = Ang[2];
3837 //delete l4a;
3838 delete l4b;
3839
3840 //StraightLine *l5a = IniPar(Vec_phi[0], Vec_Z[0], 5, Vec_phi[1], Vec_Z[1], 4);
3841 fit = Fit(line, fa - f10, 0);
3842 fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
3843 inter_5_flag = istat;
3844 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
3845 Store_Trk(fit, 5000);
3846 delete fit;
3847 StraightLine *l5b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
3848 if (Align_FLAG) {
3849 Mp->getPointAligned(2, l5b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3850 Mp->getPointAligned(3, l5b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3851 } else Mp->getPointIdealGeom(1 * 2, l5b, Up, Down, phiV_up, phiV_down);
3852 to_0_2pi(phiV_down[0]);
3853 inter_5_x = phiV_down[0];
3854 inter_5_V = phiV_down[1];
3855 inter_5_z = Down[2];
3856 HepPoint3D l1_down(cos(inter_5_x) * R_layer[1], sin(inter_5_x) * R_layer[1], inter_5_z);
3857 IncidentAngle(l5b, l1_down, pl_10->getStereoAngle(), Ang);
3858 ang_5_x = Ang[0];
3859 ang_5_v = Ang[1];
3860 ang_5 = Ang[2];
3861 //delete l5a;
3862 delete l5b;
3863 //}
3864
3865 //StraightLine *l6a = IniPar(Vec_phi[2], Vec_Z[2], 3, Vec_phi[3], Vec_Z[3], 2);
3866 fit = Fit(line, fa - f20, 0);
3867 fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
3868 inter_6_flag = istat;
3869 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
3870 Store_Trk(fit, 6000);
3871 delete fit;
3872 StraightLine *l6b = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
3873 if (Align_FLAG) {
3874 Mp->getPointAligned(4, l6b, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3875 Mp->getPointAligned(5, l6b, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3876 } else Mp->getPointIdealGeom(2 * 2, l6b, Up, Down, phiV_up, phiV_down);
3877
3878// to_0_2pi(phiV_up[0]);
3879to_0_2pi(phiV_down[0]);
3880inter_6_x = phiV_down[0];
3881inter_6_V = phiV_down[1];
3882inter_6_z = Down[2];
3883
3884HepPoint3D l2_down(cos(inter_6_x) * R_layer[2], sin(inter_6_x) * R_layer[2], inter_6_z);
3885IncidentAngle(l6b, l2_down, pl_20->getStereoAngle(), Ang);
3886ang_6_x = Ang[0];
3887ang_6_v = Ang[1];
3888ang_6 = Ang[2];
3889//delete l6a;
3890delete l6b;
3891}
3892
3893vector<double> CgemLineFit::Get_Clusters(vector<double> Vec_truth) {
3894 vector<double> clusters;
3895 StraightLine l(Vec_truth[0], Vec_truth[1], Vec_truth[2], Vec_truth[3]);
3896 double phiV_up[2], phiV_down[2];
3897 HepPoint3D Up, Down;
3898
3899 double phiV_up_Tmp[2], phiV_down_Tmp[2];
3900 HepPoint3D Up_Tmp, Down_Tmp;
3901
3902 if (!Align_Flag) Mp->getPointIdealGeom(0 * 2, l, Up, Down, phiV_up, phiV_down);
3903 else {
3904 Mp->getPointAligned(0, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3905 Mp->getPointAligned(1, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3906 }
3907 to_0_2pi(phiV_up[0]);
3908 to_0_2pi(phiV_down[0]);
3909
3910 clusters.push_back(phiV_up[0]);
3911 clusters.push_back(phiV_up[1]);
3912 clusters.push_back(phiV_down[0]);
3913 clusters.push_back(phiV_down[1]);
3914
3915 if (!Align_Flag) Mp->getPointIdealGeom(1 * 2, l, Up, Down, phiV_up, phiV_down);
3916 else {
3917 Mp->getPointAligned(2, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3918 Mp->getPointAligned(3, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3919 }
3920 to_0_2pi(phiV_up[0]);
3921 to_0_2pi(phiV_down[0]);
3922
3923 clusters.push_back(phiV_up[0]);
3924 clusters.push_back(phiV_up[1]);
3925 clusters.push_back(phiV_down[0]);
3926 clusters.push_back(phiV_down[1]);
3927
3928 if (!Align_Flag) Mp->getPointIdealGeom(2 * 2, l, Up, Down, phiV_up, phiV_down);
3929 else {
3930 Mp->getPointAligned(4, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
3931 Mp->getPointAligned(5, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
3932 }
3933 to_0_2pi(phiV_up[0]);
3934 to_0_2pi(phiV_down[0]);
3935
3936 clusters.push_back(phiV_up[0]);
3937 clusters.push_back(phiV_up[1]);
3938 clusters.push_back(phiV_down[0]);
3939 clusters.push_back(phiV_down[1]);
3940
3941 return clusters;
3942}
3943
3944/*void CgemLineFit::to_0_2pi(double &arg)
3945 {
3946 arg = fmod((arg),(2.0*TMath::Pi()));
3947 if(arg<0) arg += 2.0*TMath::Pi();
3948 else arg -= 2.0*TMath::Pi();
3949 }*/
3950void CgemLineFit::to_0_2pi(double &arg) {
3951 while (arg < -1 * acos(-1)) arg += acos(-1) * 2;
3952 while (arg > acos(-1)) arg -= acos(-1) * 2;
3953}
3954
3955double CgemLineFit::Min_diff(double arg1, double arg2) {
3956 to_0_2pi(arg1);
3957 to_0_2pi(arg2);
3958 double diff_raw = fabs(arg1 - arg2);
3959
3960 if (diff_raw > acos(-1)) return (2 * acos(-1) - diff_raw);
3961 else return diff_raw;
3962}
3963
3964double CgemLineFit::Min_dis(int R_id, double x, double v, StraightLine *l1) {
3965 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
3966 double phiV_up[2], phiV_down[2];
3967 HepPoint3D Up, Down;
3968
3969 int vir_R_id = GetVirLay(R_id, x);
3970
3971 if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
3972 else Mp->getPointAligned(vir_R_id, l1, Up, Down, phiV_up, phiV_down);
3973
3974 double ds1 = pow((phiV_up[1] - v), 2) + pow(R_layer[R_id] * Min_diff(phiV_up[0], x), 2);
3975
3976 double ds2 = pow((phiV_down[1] - v), 2) + pow(R_layer[R_id] * Min_diff(phiV_down[0], x), 2);
3977
3978 return min(ds1, ds2);
3979}
3980
3981double CgemLineFit::Min_dis_up(int R_id, double x, double z, StraightLine *l1) {
3982 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
3983 double phiV_up[2], phiV_down[2];
3984 HepPoint3D Up, Down;
3985 int vir_R_id = GetVirLay(R_id, x);
3986 if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
3987 else Mp->getPointAligned(vir_R_id, l1, Up, Down, phiV_up, phiV_down);
3988 double ds1 = pow((phiV_up[1] - z), 2) + pow(R_layer[R_id] * Min_diff(phiV_up[0], x), 2);
3989
3990 return ds1;
3991}
3992
3993double CgemLineFit::Min_dis_down(int R_id, double x, double v, StraightLine *l1) {
3994
3995 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
3996 double phiV_up[2], phiV_down[2];
3997 HepPoint3D Up, Down;
3998 int vir_R_id = GetVirLay(R_id, x);
3999 if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
4000 else Mp->getPointAligned(vir_R_id, l1, Up, Down, phiV_up, phiV_down);
4001
4002 double ds2 = pow((phiV_down[1] - v), 2) + pow(R_layer[R_id] * Min_diff(phiV_down[0], x), 2);
4003
4004 return ds2;
4005}
4006
4007bool CgemLineFit::Layer_cross(int R_id, StraightLine *l1) {
4008
4009 // need to use vir_layer instead of geo_layer -- Guoaq--2020/11/02
4010 double phiV_up[2], phiV_down[2];
4011
4012 HepPoint3D Up, Down;
4013 if (!Align_Flag) Mp->getPointIdealGeom(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
4014 else Mp->getPointAligned(R_id * 2, l1, Up, Down, phiV_up, phiV_down);
4015 if (Up[0] > 9999 || Up[1] > 9999 || Down[0] > 9999 || Down[1] > 9999) return false;
4016 else return true;
4017}
4018
4019vector<double> CgemLineFit::Trans(double par0, double par1, double par2, double par3) {
4020 vector<double> rotat;
4021 rotat.clear();
4022
4023 rotat.push_back(par0);
4024 rotat.push_back(par1);
4025 rotat.push_back(par2);
4026 rotat.push_back(par3);
4027 return rotat;
4028}
4029
4030int CgemLineFit::GetVirLay(int geolay, double phi) {
4031 int virsheet = (phi < 0) ? 0 : 1;
4032 int virlay = geolay * 2 + virsheet;
4033 return virlay;
4034}
4035
4036vector<int> CgemLineFit::GetNMaxQ(vector<double> Q_11, vector<int> L_11, int Nmax = 3) {
4037 // chose the N X clusters with the largest charge and save corresponding index
4038 vector<int> L_11_s;
4039 vector<double> Q_11_s;
4040 double t, tt;
4041 if (L_11.size() > Nmax) {
4042 for (int i = 0; i < Nmax; i++) {
4043 for (int j = i + 1; j < L_11.size(); j++) {
4044 if (Q_11[i] < Q_11[j]) {
4045 t = Q_11[i];
4046 Q_11[i] = Q_11[j];
4047 Q_11[j] = t;
4048
4049 tt = L_11[i];
4050 L_11[i] = L_11[j];
4051 L_11[j] = tt;
4052 }
4053 }
4054 Q_11_s.push_back(Q_11[i]);
4055 L_11_s.push_back(L_11[i]);
4056 }
4057
4058 } else {
4059 L_11_s.assign(L_11.begin(), L_11.end());
4060 Q_11_s.assign(Q_11.begin(), Q_11.end());
4061 }
4062 return L_11_s;
4063}
4064
4065void CgemLineFit::FilterClusters() {
4066 int max_1_0(0), max_0_0(0);
4067 int _size = Vec_flag.size();
4068 max_1_0 = _size;
4069 max_0_0 = _size;
4070
4071 for (int i = 0; i < Vec_flag.size(); i++) {
4072 if (1 == Vec_layerid[i] && 0 == Vec_sheetid[i]) {
4073 max_1_0 == i;
4074 break;
4075 }
4076 }
4077
4078 for (int j = 0; j < Vec_flag.size(); j++) {
4079 if (0 == Vec_layerid[j] && 0 == Vec_sheetid[j]) {
4080 max_0_0 == j;
4081 break;
4082 }
4083 }
4084 for (int k = Vec_flag.size() - 1; k > max_0_0; k--) { erase(k); }
4085 for (int m = Vec_flag.size() - 2; m > max_1_0; m--) { erase(m); }
4086}
4087
4088void CgemLineFit::Store_Trk(TMinuit *fit, int trkid) {
4089 double trkpar[4] = {-999, -999, -999, -999};
4090 double trkparErr[4] = {-999, -999, -999, -999};
4091
4092 Int_t ierflg, npari, nparx, istat3;
4093 Double_t fmin, edm, errdef;
4094
4095 fit->mnstat(fmin, edm, errdef, npari, nparx, istat3);
4096 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
4097 double emat[16];
4098 fit->mnemat(emat, 4);
4099
4100 RecMdcTrack *recMdcTrack = new RecMdcTrack();
4101
4102 double helixPar[5];
4103 helixPar[0] = trkpar[0];
4104 helixPar[1] = trkpar[1];
4105 helixPar[2] = 0.0;
4106 helixPar[3] = trkpar[2];
4107 helixPar[4] = trkpar[3];
4108 double errorMat[15];
4109 int k = 0;
4110 errorMat[0] = emat[0];
4111 errorMat[1] = emat[1];
4112 errorMat[2] = 0;
4113 errorMat[3] = emat[2];
4114 errorMat[4] = emat[3];
4115 errorMat[5] = emat[5];
4116 errorMat[6] = 0;
4117 errorMat[7] = emat[6];
4118 errorMat[8] = emat[7];
4119 errorMat[9] = 0;
4120 errorMat[10] = 0;
4121 errorMat[11] = 0;
4122 errorMat[12] = emat[10];
4123 errorMat[13] = emat[11];
4124 errorMat[14] = emat[15];
4125
4126 recMdcTrack->setChi2(fit->fAmin);
4127 recMdcTrack->setError(errorMat);
4128 recMdcTrack->setHelix(helixPar);
4129 recMdcTrack->setTrackId(trkid);
4130 /////
4131 //SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
4132 //RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();
4133 /////
4134 if (trkid == 0) {
4135 //cout << "before setNcluster : " << recMdcTrack->getNcluster() << endl;
4136 recMdcTrack->setNcluster(NC);
4137 //cout << "set NC = " << NC << endl;
4138 //cout << "after setNcluster : " << recMdcTrack->getNcluster() << endl;
4139 //cout << "before setVecClusters : " << recMdcTrack->getNcluster() << endl;
4140 recMdcTrack->setVecClusters(vecclusters);
4141 //cout << "after setVecClusters : " << recMdcTrack->getNcluster() << endl;
4142 m_trackList_tds->push_back(recMdcTrack);
4143 /* SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
4144 RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();
4145 cout << "after m_trackList_tds->push_back() : " << (*iter)->getNcluster() << endl; */
4146 }
4147}
4148
4149void CgemLineFit::IncidentAngle(StraightLine *cosmic_ray, HepPoint3D cluster, double theta, double ang[]) {
4150
4151 const Hep3Vector Normals(cluster[0], cluster[1], 0);
4152
4153 HepPoint3D x1 = cosmic_ray->x(1);
4154 HepPoint3D x2 = cosmic_ray->x(-1);
4155 const Hep3Vector CosmicRay(x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2]);
4156 const Hep3Vector z(0, 0, 1);
4157 const Hep3Vector phi = z.cross(Normals);
4158 Hep3Vector _phi = phi;
4159
4160 const Hep3Vector V = _phi.rotate(-1 * theta, Normals);
4161
4162 const Hep3Vector N_V = Normals.cross(V);
4163 const Hep3Vector N_phi = Normals.cross(phi);
4164
4165 Hep3Vector C_phi = CosmicRay - CosmicRay.project(N_phi);
4166 double s_C_phi = Normals.dot(C_phi);
4167 Hep3Vector v_phi = Normals.cross(C_phi);
4168 // ang[0]=Normals.angle(C_phi);
4169 if (s_C_phi <= 0) {
4170 double s_phi = v_phi.dot(z);
4171 if (s_phi >= 0)
4172 ang[0] = -(Normals.angle(-C_phi)); // 宇宙线在x-y平面上的投影和法线的夹角,空间角不分正负，所以加几个判断条件区分正负方向
4173 if (s_phi < 0) ang[0] = Normals.angle(-C_phi);
4174 /* if(s_phi==0)
4175 {
4176 ang[0]=Normals.angle(C_phi);
4177 cout<<"s_C_phi<=0 s_phi==0 "<<ang[0]<<endl;
4178 }
4179 if(ang[0]==0)
4180 cout<<"s_C_phi<=0 "<<ang[0]<<endl;*/
4181 } else if (s_C_phi > 0) {
4182 double s_phi = v_phi.dot(z);
4183 if (s_phi <= 0)
4184 ang[0] = -(Normals.angle(
4185 C_phi)); // 宇宙线在x-y平面上的投影和法线的夹角,C_phi,Normals都是向量,算出的间角不分正负，所以加几个判断条件区分正负方向
4186 if (s_phi > 0) ang[0] = Normals.angle(C_phi);
4187 /* if(s_phi==0)
4188 {
4189 ang[0]=Normals.angle(C_phi);
4190 cout<<"s_C_phi>0 s_phi==0 "<<ang[0]<<endl;
4191 }
4192 if(ang[0]==0)
4193 cout<<"s_C_phi>0 "<<ang[0]<<endl;*/
4194 }
4195 Hep3Vector C_V = CosmicRay - CosmicRay.project(N_V); // C_V是宇宙线在与V向垂直的平面上的投影
4196 double s_C_V = Normals.dot(C_V);
4197 Hep3Vector v_V = Normals.cross(C_V);
4198 // ang[1]=Normals.angle(C_V);
4199 if (s_C_V <= 0) {
4200 double s_V = v_V.dot(z);
4201 if (s_V >= 0) ang[1] = -(Normals.angle(-C_V)); // 宇宙线在与V向垂直的平面上的投影和法线的夹角
4202 if (s_V < 0) ang[1] = Normals.angle(-C_V);
4203 /* if(s_V==0)
4204 {
4205 ang[1]=Normals.angle(C_V);
4206 cout<<"s_C_V<=0 s_V==0 "<<ang[1]<<endl;
4207 }
4208 if(ang[1]==0)
4209 cout<<"s_C_V<=0 "<<ang[1]<<endl;*/
4210 } else if (s_C_V > 0) {
4211 double s_V = v_V.dot(z);
4212 if (s_V <= 0) ang[1] = -(Normals.angle(C_V));
4213 if (s_V > 0) ang[1] = Normals.angle(C_V);
4214 /* if(s_V==0)
4215 {
4216 ang[1]=Normals.angle(C_V);
4217 cout<<"s_C_V>0 s_V==0 "<<ang[1]<<endl;
4218 }
4219 if(ang[1]==0)
4220 cout<<"s_C_V>0 "<<ang[1]<<endl;*/
4221 }
4222 ang[2] = Normals.angle(CosmicRay); // 宇宙线和法线的夹角即是入射角,但这是个空间上的三维的角，所以我们要用投影平面的夹角。
4223 /*Hep3Vector C_V=CosmicRay-CosmicRay.project(N_V);
4224 ang[1]=Normals.angle(C_V);
4225 ang[2]=Normals.angle(CosmicRay);*/
4226}
4227
4228void CgemLineFit::GetRes(StraightLine *l, int i_inter) {
4229 int i_layer = 0;
4230 int i_sheet = 0;
4231 if (i_inter == 1 || i_inter == 6) {
4232 i_layer = 2;
4233 if (i_inter == 1) i_sheet = 1;
4234 }
4235 if (i_inter == 2 || i_inter == 5) {
4236 i_layer = 1;
4237 if (i_inter == 2) i_sheet = 1;
4238 }
4239 CgemGeoReadoutPlane *readoutPlane;
4240 readoutPlane = m_cgemGeomSvc->getReadoutPlane(i_layer, i_sheet);
4241
4242 double phiV_up[2], phiV_down[2];
4243 HepPoint3D Up, Down;
4244 double phiV_up_Tmp[2], phiV_down_Tmp[2];
4245 HepPoint3D Up_Tmp, Down_Tmp;
4246
4247 if (!Align_Flag) Mp->getPointIdealGeom(i_layer * 2, l, Up, Down, phiV_up, phiV_down);
4248 else {
4249 Mp->getPointAligned(i_layer * 2, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
4250 Mp->getPointAligned(i_layer * 2 + 1, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
4251 }
4252 double V, phi, min_dst(99999), min_V(99999), min_phi(99999);
4253 double Z, min_Z(99999);
4254
4255 V = phiV_up[1];
4256 phi = phiV_up[0];
4257 //Z = Up[2];
4258
4259 if (i_inter == 4 || i_inter == 5 || i_inter == 6) {
4260 V = phiV_down[1];
4261 phi = phiV_down[0];
4262 //Z = Down[2];
4263 }
4264 Z=readoutPlane->getZFromPhiV(phi, V);
4265
4266 int _loop = 0;
4267 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
4268 RecCgemClusterCol::iterator iter = recCgemClusterCol->begin();
4269 int ind = 0;
4270 for (; iter != recCgemClusterCol->end(); iter++) {
4271 _loop++;
4272 RecCgemCluster *cluster = (*iter);
4273 int flag = cluster->getflag();
4274
4275 int layerid = cluster->getlayerid();
4276 if (layerid != i_layer) continue;
4277 int sheetid = cluster->getsheetid();
4278 if (sheetid != i_sheet) continue;
4279 double _Q = cluster->getenergydeposit();
4280 double _phi = cluster->getrecphi();
4281 double _V = cluster->getrecv();
4282 double _Z = cluster->getRecZ();
4283 double _size = fabs(cluster->getclusterflagb() - cluster->getclusterflage());
4284 if (flag != 2) continue;
4285 double dV = _V - V;
4286 double dX = (_phi - phi) * R_layer[layerid];
4287 double dst = pow(dV, 2.0) + pow(dX, 2.0) + 2 * abs(dX * dV) * cos(M_PI * AngleOfStereo[layerid] / 180);
4288 // if(event==6)
4289 // cout<<__LINE__<<" dV "<<dV<<" dX "<<dX<<" dst "<<dst<<" _phi
4290 // "<<_phi<<" _V "<<_V<<" phi "<<phi<<" V "<<V<<endl;
4291 if (abs(dst) < min_dst) {
4292 min_dst = dst;
4293 min_V = _V;
4294 min_phi = _phi;
4295 min_Z = _Z;
4296 }
4297 }
4298
4299 Test_phi = phi;
4300 Test_V = V;
4301 Test_Z = Z;
4302 Rec_phi = min_phi;
4303 Rec_V = min_V;
4304 Rec_Z = min_Z;
4305
4306
4307 double Ang[3];
4308 if (i_inter > 0) {
4309
4310 Double_t fmin, edm, errdef;
4311 Int_t ierflg, npari, nparx, istat;
4312 double trkpar[4] = {-999, -999, -999, -999};
4313 double trkparErr[4] = {-999, -999, -999, -999};
4314 int test_flag, OC; // orderclusters
4315 double StereoAngle;
4316
4317 // 0 1 2 3 4 5
4318 // order is layer 2 , top, down, layer 1, top down, layer 0, top down
4319
4320 switch (i_inter) {
4321 case 1:
4322 test_flag = f21;
4323 OC = 0;
4324 StereoAngle = pl_21->getStereoAngle();
4325 break;
4326 case 2:
4327 test_flag = f11;
4328 OC = 2;
4329 StereoAngle = pl_11->getStereoAngle();
4330 break;
4331 case 3:
4332 test_flag = f01;
4333 OC = 4;
4334 StereoAngle = pl_00->getStereoAngle();
4335 break;
4336 case 4:
4337 test_flag = f00;
4338 OC = 5;
4339 StereoAngle = pl_00->getStereoAngle();
4340 break;
4341 case 5:
4342 test_flag = f10;
4343 OC = 3;
4344 StereoAngle = pl_10->getStereoAngle();
4345 break;
4346 case 6:
4347 test_flag = f20;
4348 OC = 1;
4349 StereoAngle = pl_20->getStereoAngle();
4350 break;
4351 default: break;
4352 }
4353
4354 TMinuit *fit = Fit(l, fa - test_flag, 0);
4355 fit->mnstat(fmin, edm, errdef, npari, nparx, istat);
4356 inter_flag = istat;
4357 for (int i = 0; i < 4; i++) { fit->GetParameter(i, trkpar[i], trkparErr[i]); }
4358 StraightLine *aline = new StraightLine(trkpar[0], trkpar[1], trkpar[2], trkpar[3]);
4359 // Store_Trk(fit, 1000);
4360 delete fit;
4361
4362 if (!Align_Flag) Mp->getPointIdealGeom(i_layer * 2, l, Up, Down, phiV_up, phiV_down);
4363 else {
4364 Mp->getPointAligned(i_layer * 2, l, Up_Tmp, Down, phiV_up_Tmp, phiV_down);
4365 Mp->getPointAligned(i_layer * 2 + 1, l, Up, Down_Tmp, phiV_up, phiV_down_Tmp);
4366 }
4367 to_0_2pi(phiV_up[0]);
4368 to_0_2pi(phiV_down[0]);
4369
4370 inter_x = phiV_up[0];
4371 inter_V = phiV_up[1];
4372 inter_z = Up[2];
4373
4374 if (i_inter == 4 || i_inter == 5 || i_inter == 6) {
4375 inter_x = phiV_down[0];
4376 inter_V = phiV_down[1];
4377 inter_z = Down[2];
4378 }
4379
4380 //cout << "OC = " << OC << endl;
4381
4382 //HepPoint3D point(cos(Vec_phi[OC]) * R_layer[i_layer], sin(Vec_phi[OC]) * R_layer[i_layer], Vec_Z[OC]);
4383 //cout << "phi = " << Vec_phi[OC] << endl;
4384 HepPoint3D point(cos(inter_x) * R_layer[i_layer], sin(inter_x) * R_layer[i_layer], inter_z);
4385 IncidentAngle(aline, point, StereoAngle, Ang);
4386 ang_x = Ang[0];
4387 ang_v = Ang[1];
4388 Angle = Ang[2];
4389 //cout << "inter_x = " << inter_x << " ang_x = " << ang_x << endl;
4390
4391 }
4392}
4393
4394double CgemLineFit::RealPhi(double SimuPhi) {
4395 double RPhi;
4396 if (SimuPhi <= TMath::Pi()) RPhi = SimuPhi;
4397 else RPhi = SimuPhi - 2 * TMath::Pi();
4398 return RPhi;
4399}
4400
4401double CgemLineFit::CalSigma2(double layerid, double flag, double x, double V, double z, StraightLine *cosmic_ray) {
4402 double sigma = 0;
4403 double Ang[3];
4404 if (layerid == 0) {
4405 HepPoint3D l0(cos(x) * R_layer[0], sin(x) * R_layer[0], z);
4406 IncidentAngle(cosmic_ray, l0, pl_00->getStereoAngle(), Ang);
4407 if (flag == 0) {
4408 const int n = 28;
4409 double angle[n] = {-1.35, -1.25, -1.15, -1.05, -0.95, -0.85, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05,
4410 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35};
4411 double y[n] = {2.3576, 1.44448, 1.137, 0.998306, 0.819508, 0.695739, 0.583851, 0.487318, 0.408487, 0.317606,
4412 0.245717, 0.178443, 0.109751, 0.0469862, 0.0463951, 0.109757, 0.180819, 0.252827, 0.326331, 0.409922,
4413 0.482766, 0.594053, 0.690917, 0.835596, 0.93621, 1.25398, 1.56933, 1.93546};
4414 double *ex = nullptr;
4415 double ey[n] = {0.179035, 0.0641824, 0.0347424, 0.0256942, 0.015515, 0.0123214, 0.00893213,
4416 0.00584735, 0.00470714, 0.00334024, 0.00247131, 0.00171306, 0.00101272, 0.000571988,
4417 0.000525056, 0.00106782, 0.00170087, 0.00260604, 0.00347387, 0.00486819, 0.00638214,
4418 0.00910626, 0.0117511, 0.0172709, 0.0236933, 0.040096, 0.0607565, 0.171064};
4419 TGraphErrors gr(n, angle, y, ex, ey);
4420 sigma = gr.Eval(Ang[0]);
4421 }
4422 if (flag == 1) {
4423 const int n = 26;
4424 double angle[n] = {-1.25, -1.15, -1.05, -0.95, -0.85, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05,
4425 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25};
4426 double y[n] = {1.32155, 1.00793, 0.844239, 0.745985, 0.607887, 0.525712, 0.449061, 0.372476, 0.305783,
4427 0.231835, 0.167037, 0.105903, 0.0526632, 0.0518158, 0.106611, 0.170067, 0.235356, 0.301801,
4428 0.374196, 0.46125, 0.526767, 0.61776, 0.749122, 0.848473, 1.01892, 1.50534};
4429 double *ex = nullptr;
4430 double ey[n] = {0.108111, 0.0523342, 0.0275105, 0.0227304, 0.0135576, 0.0095431, 0.00675217,
4431 0.00439062, 0.00328685, 0.00225977, 0.0013645, 0.000902477, 0.000515666, 0.000505206,
4432 0.000917112, 0.00148768, 0.00221462, 0.0031313, 0.0046925, 0.00745832, 0.00983847,
4433 0.0151778, 0.0227792, 0.033388, 0.0477401, 0.159787};
4434 TGraphErrors gr(n, angle, y, ex, ey);
4435 sigma = gr.Eval(Ang[1]);
4436 }
4437 }
4438 if (layerid == 1) {
4439 HepPoint3D l1(cos(x) * R_layer[1], sin(x) * R_layer[1], z);
4440 IncidentAngle(cosmic_ray, l1, pl_11->getStereoAngle(), Ang);
4441 if (flag == 0) {
4442 const int n = 14;
4443 double angle[n] = {-0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65};
4444 double y[n] = {0.459771, 0.381421, 0.310736, 0.238213, 0.169008, 0.103486, 0.0461908,
4445 0.0453169, 0.103585, 0.166698, 0.240332, 0.3125, 0.37973, 0.462382};
4446 double *ex = nullptr;
4447 double ey[n] = {0.00685022, 0.00437997, 0.00320362, 0.00214694, 0.00136347, 0.00081222, 0.000440626,
4448 0.000435363, 0.000818715, 0.00128686, 0.00214114, 0.00313311, 0.00436737, 0.0071195};
4449
4450 TGraphErrors gr(n, angle, y, ex, ey);
4451 sigma = gr.Eval(Ang[0]);
4452 }
4453 if (flag == 1) {
4454 const int n = 16;
4455 double angle[n] = {-0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75};
4456 double y[n] = {0.543351, 0.449597, 0.376926, 0.299783, 0.234983, 0.165406, 0.105945, 0.0544188,
4457 0.0540844, 0.10723, 0.167146, 0.234566, 0.300303, 0.37238, 0.461751, 0.570197};
4458 double *ex = nullptr;
4459 double ey[n] = {0.0255211, 0.00933694, 0.00529557, 0.00317545, 0.00209458, 0.0013353, 0.000773631, 0.00046874,
4460 0.000453118, 0.000828085, 0.00125195, 0.00202037, 0.00315857, 0.00489179, 0.0100004, 0.032045};
4461 TGraphErrors gr(n, angle, y, ex, ey);
4462 sigma = gr.Eval(Ang[1]);
4463 }
4464 }
4465 /////
4466 if (layerid == 2) {
4467 HepPoint3D l2(cos(x) * R_layer[2], sin(x) * R_layer[2], z);
4468 IncidentAngle(cosmic_ray, l2, pl_21->getStereoAngle(), Ang);
4469 if (flag == 0) {
4470 const int n = 14;
4471 double angle[n] = {-0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65};
4472 double y[n] = {0.459771, 0.381421, 0.310736, 0.238213, 0.169008, 0.103486, 0.0461908,
4473 0.0453169, 0.103585, 0.166698, 0.240332, 0.3125, 0.37973, 0.462382};
4474 double *ex = nullptr;
4475 double ey[n] = {0.00685022, 0.00437997, 0.00320362, 0.00214694, 0.00136347, 0.00081222, 0.000440626,
4476 0.000435363, 0.000818715, 0.00128686, 0.00214114, 0.00313311, 0.00436737, 0.0071195};
4477
4478 TGraphErrors gr(n, angle, y, ex, ey);
4479 sigma = gr.Eval(Ang[0]);
4480 }
4481 if (flag == 1) {
4482 const int n = 16;
4483 double angle[n] = {-0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75};
4484 double y[n] = {0.543351, 0.449597, 0.376926, 0.299783, 0.234983, 0.165406, 0.105945, 0.0544188,
4485 0.0540844, 0.10723, 0.167146, 0.234566, 0.300303, 0.37238, 0.461751, 0.570197};
4486 double *ex = nullptr;
4487 double ey[n] = {0.0255211, 0.00933694, 0.00529557, 0.00317545, 0.00209458, 0.0013353, 0.000773631, 0.00046874,
4488 0.000453118, 0.000828085, 0.00125195, 0.00202037, 0.00315857, 0.00489179, 0.0100004, 0.032045};
4489 TGraphErrors gr(n, angle, y, ex, ey);
4490 sigma = gr.Eval(Ang[1]);
4491 }
4492 }
4493 /////
4494 return sigma * sigma;
4495}
4496double CgemLineFit::CalLineSAtY(StraightLine &l, double y) {
4497 double drho = l.dr();
4498 double phi0 = l.phi0();
4499 double s = (y - drho * sin(phi0)) / sin(phi0 - TMath::Pi() / 2);
4500 return s;
4501}
4502float CgemLineFit::get_Time(CgemDigiCol::iterator iDigiCol) {
4503 // Get digi time
4504 float time = (*iDigiCol)->getTime_ns();
4505 // Get rising time from calibration
4506 float time_rising = get_TimeRising(iDigiCol);
4507 // Get time-walk
4508 float time_walk = get_TimeWalk(iDigiCol);
4509 // Get time-reference
4510 float time_reference = get_TimeReference(iDigiCol);
4511 //
4512 float time_shift_custom = -35;
4513 // cout<<time<<" "<<time_rising<<" "<<time_walk<<" "<<time_reference;
4514 time -= (time_rising + time_walk + time_reference + time_shift_custom);
4515 // cout<<" "<<time<<endl;
4516 return time;
4517}
4518
4519float CgemLineFit::get_TimeRising(CgemDigiCol::iterator iDigiCol) {
4520 float time_rising = 0;
4521 // Get the digi information
4522 const Identifier ident = (*iDigiCol)->identify();
4523 int layerid = CgemID::layer(ident);
4524 int sheetid = CgemID::sheet(ident);
4525 int stripid = CgemID::strip(ident);
4526 int view = 1;
4527 bool is_xstrip = CgemID::is_xstrip(ident);
4528 if (is_xstrip == true) view = 0;
4529 float charge = (*iDigiCol)->getCharge_fc();
4530 // Get rising time from calibration
4531 time_rising = myCgemCalibSvc->getTimeRising(layerid, view, sheetid, stripid, charge, 0.);
4532 return time_rising;
4533}
4534
4535float CgemLineFit::get_TimeWalk(CgemDigiCol::iterator iDigiCol) {
4536 float time_walk = 0;
4537 const Identifier ident = (*iDigiCol)->identify();
4538 int layerid = CgemID::layer(ident);
4539 int sheetid = CgemID::sheet(ident);
4540 int stripid = CgemID::strip(ident);
4541 int view = 1;
4542 bool is_xstrip = CgemID::is_xstrip(ident);
4543 if (is_xstrip == true) view = 0;
4544 float thr = lutreader->Get_thr_T_fC(layerid, sheetid, view, stripid);
4545 float charge = (*iDigiCol)->getChargeChannel();
4546 time_walk = myCgemCalibSvc->getTimeWalk(charge, thr);
4547 return time_walk;
4548}
4549
4550float CgemLineFit::get_TimeReference(CgemDigiCol::iterator iDigiCol) {
4551 float time_reference = 0;
4552 const Identifier ident = (*iDigiCol)->identify();
4553 int layerid = CgemID::layer(ident);
4554 int sheetid = CgemID::sheet(ident);
4555 int stripid = CgemID::strip(ident);
4556 int view = 1;
4557 bool is_xstrip = CgemID::is_xstrip(ident);
4558 if (is_xstrip == true) view = 0;
4559 time_reference += lutreader->GetSignal_FEBStartTime_ns(layerid, sheetid, view, stripid);
4560 time_reference += lutreader->GetSignal_StartTime_ns(layerid, sheetid, view, stripid);
4561 return time_reference;
4562}
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:100
Vec_m_layer
vector< double > Vec_m_layer
Definition CgemLineFit.cxx:97
Vec_Virlayerid
vector< double > Vec_Virlayerid
Definition CgemLineFit.cxx:91
Vec_VstripQ
vector< vector< double > > Vec_VstripQ
Definition CgemLineFit.cxx:94
f21
int f21(1)
_clst_2
int _clst_2(0)
phi0_set
double phi0_set
Definition CgemLineFit.cxx:121
GeoLayer2
CgemGeoLayer * GeoLayer2
Definition CgemLineFit.cxx:103
Vec_Q
vector< double > Vec_Q
Definition CgemLineFit.cxx:91
Vec_x_stripsid
vector< int > Vec_x_stripsid
Definition CgemLineFit.cxx:95
f11
int f11(2)
Vec_flag
vector< double > Vec_flag
Definition CgemLineFit.cxx:91
sigma2
double sigma2(0)
f01
int f01(8)
Vec_VstripT
vector< vector< double > > Vec_VstripT
Definition CgemLineFit.cxx:94
AngleOfStereo
double AngleOfStereo[3]
Definition CgemLineFit.cxx:123
R_layer
double R_layer[3]
Definition CgemLineFit.cxx:122
Vec_x_stripsTf
vector< double > Vec_x_stripsTf
Definition CgemLineFit.cxx:93
Vec_m_layerid
vector< double > Vec_m_layerid
Definition CgemLineFit.cxx:97
Vec_XstripT
vector< vector< double > > Vec_XstripT
Definition CgemLineFit.cxx:94
Vec_Z
vector< double > Vec_Z
Definition CgemLineFit.cxx:91
f20
int f20(32)
Vec_m_phi
vector< double > Vec_m_phi
Definition CgemLineFit.cxx:97
Vec_XQ
vector< double > Vec_XQ
Definition CgemLineFit.cxx:92
Vec_VstripTf
vector< vector< double > > Vec_VstripTf
Definition CgemLineFit.cxx:94
Vec_VstripID
vector< vector< int > > Vec_VstripID
Definition CgemLineFit.cxx:96
Mp
CgemMidDriftPlane * Mp
Definition CgemLineFit.cxx:99
Vec_XstripTf
vector< vector< double > > Vec_XstripTf
Definition CgemLineFit.cxx:94
Vec_v_stripsQ
vector< double > Vec_v_stripsQ
Definition CgemLineFit.cxx:93
Vec_x_stripsT
vector< double > Vec_x_stripsT
Definition CgemLineFit.cxx:93
fa
int fa(63)
Vec_m_Q
vector< double > Vec_m_Q
Definition CgemLineFit.cxx:97
Vec_m_Z
vector< double > Vec_m_Z
Definition CgemLineFit.cxx:97
f00
int f00(4)
length
double length
Definition CgemLineFit.cxx:121
GeoLayer0
CgemGeoLayer * GeoLayer0
Definition CgemLineFit.cxx:101
Vec_sheetid
vector< double > Vec_sheetid
Definition CgemLineFit.cxx:91
f10
int f10(16)
Vec_m_v
vector< double > Vec_m_v
Definition CgemLineFit.cxx:97
Vec_v_stripsTf
vector< double > Vec_v_stripsTf
Definition CgemLineFit.cxx:93
l_outer
StraightLine * l_outer
Definition CgemLineFit.cxx:104
pl_00
CgemGeoReadoutPlane * pl_00
Definition CgemLineFit.cxx:106
_clst_0
int _clst_0(0)
Vec_xCluSize
vector< double > Vec_xCluSize
Definition CgemLineFit.cxx:91
Vec_VQ
vector< double > Vec_VQ
Definition CgemLineFit.cxx:92
Vec_v_stripsid
vector< int > Vec_v_stripsid
Definition CgemLineFit.cxx:95
tanl_set
double tanl_set
Definition CgemLineFit.cxx:121
Vec_vCluSize
vector< double > Vec_vCluSize
Definition CgemLineFit.cxx:92
NC
int NC
Definition CgemLineFit.cxx:124
GeoLayer1
CgemGeoLayer * GeoLayer1
Definition CgemLineFit.cxx:102
Vec_v_stripsT
vector< double > Vec_v_stripsT
Definition CgemLineFit.cxx:93
Vec_layer
vector< double > Vec_layer
Definition CgemLineFit.cxx:91
Vec_m_flag
vector< double > Vec_m_flag
Definition CgemLineFit.cxx:97
Vec_v
vector< double > Vec_v
Definition CgemLineFit.cxx:91
Vec_XstripQ
vector< vector< double > > Vec_XstripQ
Definition CgemLineFit.cxx:94
_clst_1
int _clst_1(0)
pl_20
CgemGeoReadoutPlane * pl_20
Definition CgemLineFit.cxx:109
dz_set
double dz_set
Definition CgemLineFit.cxx:121
Align_FLAG
bool Align_FLAG(false)
Vec_phi
vector< double > Vec_phi
Definition CgemLineFit.cxx:91
sheet_flag
int sheet_flag(0)
Vec_layerid
vector< double > Vec_layerid
Definition CgemLineFit.cxx:91
Vec_x_stripsQ
vector< double > Vec_x_stripsQ
Definition CgemLineFit.cxx:93
Vec_XstripID
vector< vector< int > > Vec_XstripID
Definition CgemLineFit.cxx:96
dr_set
double dr_set
Definition CgemLineFit.cxx:121
pl_11
CgemGeoReadoutPlane * pl_11
Definition CgemLineFit.cxx:108
pl_10
CgemGeoReadoutPlane * pl_10
Definition CgemLineFit.cxx:107
pl_21
CgemGeoReadoutPlane * pl_21
Definition CgemLineFit.cxx:110
Vec_m_sheetid
vector< double > Vec_m_sheetid
Definition CgemLineFit.cxx:97
n
const Int_t n
Definition DataBase/tau_mode.c:65
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
gr
TGraph * gr
Definition DataBase/tau_mode.c:69
time
Double_t time
Definition DataBase/tau_mode.c:7
arg
double arg(const EvtComplex &c)
Definition EvtComplex.hh:227
abs
double abs(const EvtComplex &c)
Definition EvtComplex.hh:212
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
s
XmlRpcServer s
Definition HelloServer.cpp:11
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
M_PI
#define M_PI
Definition TConstant.h:4
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::getInnerROfAnodeCu1
double getInnerROfAnodeCu1() const
Definition CgemGeoLayer.h:163
CgemGeoLayer::getNumberOfSheet
int getNumberOfSheet() const
Definition CgemGeoLayer.h:215
CgemGeoLayer::getInnerROfAnodeCu2
double getInnerROfAnodeCu2() const
Definition CgemGeoLayer.h:167
CgemGeoReadoutPlane::getVFromPhiZ
double getVFromPhiZ(double phi, double z, bool checkRange=true) const
Definition CgemGeoReadoutPlane.h:178
CgemGeoReadoutPlane::getZFromPhiV
double getZFromPhiV(double phi, double V, int checkXRange=1) const
Definition CgemGeoReadoutPlane.h:213
CgemGeoReadoutPlane::OnThePlane
bool OnThePlane(double phi, double z) const
Definition CgemGeoReadoutPlane.h:154
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition CgemGeomSvc.h:48
CgemGeomSvc::getLengthOfCgem
double getLengthOfCgem() const
Definition CgemGeomSvc.h:42
CgemGeomSvc::getNumberOfCgemLayer
int getNumberOfCgemLayer() const
Definition CgemGeomSvc.h:45
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
CgemID::strip
static int strip(const Identifier &id)
Definition CgemID.cxx:83
CgemID::sheet
static int sheet(const Identifier &id)
Definition CgemID.cxx:77
CgemID::layer
static int layer(const Identifier &id)
Definition CgemID.cxx:71
CgemID::is_xstrip
static bool is_xstrip(const Identifier &id)
Definition CgemID.cxx:64
CgemLUTReader::GetSignal_StartTime_ns
float GetSignal_StartTime_ns(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:342
CgemLUTReader::Get_thr_T_fC
float Get_thr_T_fC(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:282
CgemLUTReader::GetSignal_FEBStartTime_ns
float GetSignal_FEBStartTime_ns(int ilayer, int isheet, int iview, int istrip)
Definition CgemLUTReader.cxx:360
CgemLineFit::Rec_Clusters_mTPC
void Rec_Clusters_mTPC()
Definition CgemLineFit.cxx:3462
CgemLineFit::GetNMaxQ
static vector< int > GetNMaxQ(vector< double > Q_11, vector< int > L_11, int Nmax)
Definition CgemLineFit.cxx:4036
CgemLineFit::finalize
StatusCode finalize()
Definition CgemLineFit.cxx:785
CgemLineFit::OrderClusters
void OrderClusters()
Definition CgemLineFit.cxx:2834
CgemLineFit::Rec_ClusterQ
void Rec_ClusterQ()
Definition CgemLineFit.cxx:3528
CgemLineFit::OrderClusterSizeQ
void OrderClusterSizeQ()
Definition CgemLineFit.cxx:2880
CgemLineFit::Filter
void Filter(int layerid, StraightLine *l1)
Definition CgemLineFit.cxx:3008
CgemLineFit::erase
void erase(int i)
Definition CgemLineFit.cxx:2955
CgemLineFit::CalSigma2
static double CalSigma2(double layerid, double flag, double x, double V, double z, StraightLine *cosmic_ray)
Definition CgemLineFit.cxx:4401
CgemLineFit::Swap
void Swap(int i, int j)
Definition CgemLineFit.cxx:2863
CgemLineFit::IniPar
StraightLine * IniPar(double phi1, double z1, int i, double phi2, double z2, int j)
Definition CgemLineFit.cxx:2991
CgemLineFit::Get4Par
vector< double > Get4Par(HepLorentzVector p4, Hep3Vector bp)
Definition CgemLineFit.cxx:3319
CgemLineFit::Min_dis_up
double Min_dis_up(int R_id, double x, double z, StraightLine *l1)
Definition CgemLineFit.cxx:3981
CgemLineFit::execute
StatusCode execute()
Definition CgemLineFit.cxx:561
CgemLineFit::CgemLineFit
CgemLineFit(const std::string &name, ISvcLocator *pSvcLocator)
Definition CgemLineFit.cxx:127
CgemLineFit::GetVirLay
int GetVirLay(int geolay, double phi)
Definition CgemLineFit.cxx:4030
CgemLineFit::Fit_Clusters
void Fit_Clusters(double par[])
Definition CgemLineFit.cxx:3680
CgemLineFit::fcn2
static void fcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition CgemLineFit.cxx:2644
CgemLineFit::GetIntersection
void GetIntersection(StraightLine *line)
Definition CgemLineFit.cxx:3705
CgemLineFit::Layer_cross
bool Layer_cross(int R_id, StraightLine *l1)
Definition CgemLineFit.cxx:4007
CgemLineFit::Trans
static vector< double > Trans(double par0, double par1, double par2, double par3)
Definition CgemLineFit.cxx:4019
CgemLineFit::Min_dis_down
double Min_dis_down(int R_id, double x, double z, StraightLine *l1)
Definition CgemLineFit.cxx:3993
CgemLineFit::fcn3
static void fcn3(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition CgemLineFit.cxx:2668
CgemLineFit::fcn
static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition CgemLineFit.cxx:2623
CgemLineFit::IncidentAngle
static void IncidentAngle(StraightLine *cosmic_ray, HepPoint3D cluster, double theta, double ang[])
Definition CgemLineFit.cxx:4149
CgemLineFit::dV
static double dV(int layer, double R, double dr, double phi0, double z0, double tglam, double x, double V)
Definition CgemLineFit.cxx:2720
CgemLineFit::Min_dis
double Min_dis(int R_id, double x, double z, StraightLine *l1)
Definition CgemLineFit.cxx:3964
CgemLineFit::MC_truth
vector< double > MC_truth()
Definition CgemLineFit.cxx:3136
CgemLineFit::HitPos
bool HitPos(HepLorentzVector p4, Hep3Vector bp)
Definition CgemLineFit.cxx:3342
CgemLineFit::Switch_CCmTPC
int Switch_CCmTPC
Definition CgemLineFit.h:68
CgemLineFit::to_0_2pi
static void to_0_2pi(double &arg)
Definition CgemLineFit.cxx:3950
CgemLineFit::Data_Closest
bool Data_Closest()
Definition CgemLineFit.cxx:1180
CgemLineFit::Fit
TMinuit * Fit(StraightLine *l_outer, int sheetflag, int typ)
Definition CgemLineFit.cxx:2738
CgemLineFit::dx
static double dx(double layerid, double R, double dr, double phi0, double z0, double tanl, double x)
Definition CgemLineFit.cxx:2695
CgemLineFit::Get_OtherIniPar
void Get_OtherIniPar(int clst_0, int clst_1, int clst_2)
Definition CgemLineFit.cxx:3367
CgemLineFit::FilterClusters
void FilterClusters()
Definition CgemLineFit.cxx:4065
CgemLineFit::Min_diff
static double Min_diff(double arg1, double arg2)
Definition CgemLineFit.cxx:3955
CgemLineFit::initialize
StatusCode initialize()
Definition CgemLineFit.cxx:147
CgemLineFit::Store_Trk
void Store_Trk(TMinuit *fit, int trkid)
Definition CgemLineFit.cxx:4088
CgemLineFit::CalLineSAtY
static double CalLineSAtY(StraightLine &l, double y)
Definition CgemLineFit.cxx:4496
CgemLineFit::Discard
void Discard(int layer)
Definition CgemLineFit.cxx:2985
CgemLineFit::Rec_ClusterSize
void Rec_ClusterSize()
Definition CgemLineFit.cxx:3499
CgemLineFit::MinQ_Clus2D
double MinQ_Clus2D
Definition CgemLineFit.h:120
CgemLineFit::RealPhi
double RealPhi(double SimuPhi)
Definition CgemLineFit.cxx:4394
CgemLineFit::Get_Clusters
vector< double > Get_Clusters(vector< double >Vec_truth)
Definition CgemLineFit.cxx:3893
CgemLineFit::Rec_Clusters
void Rec_Clusters()
Definition CgemLineFit.cxx:3422
CgemLineFit::registerTrack
StatusCode registerTrack(RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
Definition CgemLineFit.cxx:2574
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
ICgemCalibFunSvc::getTimeRising
virtual double getTimeRising(int layer, int xvFlag, int sheet, int stripID, double Q=100., double z=0.) const =0
ICgemCalibFunSvc::getTimeWalk
virtual double getTimeWalk(int layer, int xvFlag, int sheet, int stripID, double Q) const =0
RecCgemCluster::getrecphi_mTPC
double getrecphi_mTPC(void) const
Definition RecCgemCluster.h:66
RecCgemCluster::getenergydeposit
double getenergydeposit(void) const
Definition RecCgemCluster.h:59
RecCgemCluster::getRecZ
double getRecZ(void) const
Definition RecCgemCluster.h:62
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition RecCgemCluster.h:54
RecCgemCluster::getrecv_mTPC
double getrecv_mTPC(void) const
Definition RecCgemCluster.h:67
RecCgemCluster::getlayerid
int getlayerid(void) const
Definition RecCgemCluster.h:56
RecCgemCluster::getflag
int getflag(void) const
Definition RecCgemCluster.h:58
RecCgemCluster::getclusterflagb
int getclusterflagb(void) const
Definition RecCgemCluster.h:73
RecCgemCluster::getRecZ_mTPC
double getRecZ_mTPC(void) const
Definition RecCgemCluster.h:68
RecCgemCluster::getrecphi
double getrecphi(void) const
Definition RecCgemCluster.h:60
RecCgemCluster::getrecv
double getrecv(void) const
Definition RecCgemCluster.h:61
RecCgemCluster::setTrkId
void setTrkId(int trkid)
Definition RecCgemCluster.h:36
RecCgemCluster::getsheetid
int getsheetid(void) const
Definition RecCgemCluster.h:57
RecCgemCluster::getclusterflage
int getclusterflage(void) const
Definition RecCgemCluster.h:74
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