CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
XtInteCalib.cpp
Go to the documentation of this file.
1#include "include/XtInteCalib.h"
2#include "include/fun.h"
3
4#include <cmath>
5
6#include "TMinuit.h"
7#include "TF1.h"
8#include "TTree.h"
9
10XtInteCalib::XtInteCalib(){
11 m_nMaxGrPoint = 5000;
12 for(int lay=0; lay<NLAYER; lay++){
13 m_tbinWid[lay][0] = 5.0;
14 m_tbinWid[lay][1] = 10.0;
15 m_tbinWid[lay][2] = 20.0;
16
17 m_tbinLim[lay][0] = -10.0;
18 m_tbinLim[lay][1] = 30.0;
19 if(lay < 8) m_tbinLim[lay][2] = 210.0;
20 else m_tbinLim[lay][2] = 350.0;
21 m_tbinLim[lay][3] = 990.0;
22 }
23 cout << "Calibration type: XtInteCalib" << endl;
24}
25
28
29void XtInteCalib::init(TObjArray* hlist, MdcCosGeom* pGeom){
30 CalibBase::init(hlist, pGeom);
31
32 m_fdPf = new TFolder("mfdProfile", "fdProfile");
33 hlist -> Add(m_fdPf);
34
35 m_haxis = new TH2F("axis", "", 200, -50, 1000, 50, 0, 10);
36 m_haxis -> SetStats(0);
37 m_fdPf -> Add(m_haxis);
38
39 char hname[200];
40 for(int lay=0; lay<NLAYER; lay++){
41 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
42 for(int lr=0; lr<2; lr++){
43 sprintf(hname, "mxt%02d_%02d_%d_gr", lay, iEntr, lr);
44 m_grXt[lay][iEntr][lr] = new TGraph();
45 m_grXt[lay][iEntr][lr]->SetName(hname);
46 m_grXt[lay][iEntr][lr]->SetMarkerColor(2);
47 m_fdPf->Add(m_grXt[lay][iEntr][lr]);
48
49 int xbinN = (int)((m_tbinLim[lay][1] - m_tbinLim[lay][0])/m_tbinWid[lay][0] + 0.5);
50 sprintf(hname, "mxt%02d_%02d_%d_near", lay, iEntr, lr);
51 m_pfNear[lay][iEntr][lr] = new TProfile(hname, hname, xbinN, m_tbinLim[lay][0], m_tbinLim[lay][1]);
52 m_fdPf->Add(m_pfNear[lay][iEntr][lr]);
53
54 int xbinM = (int)((m_tbinLim[lay][2] - m_tbinLim[lay][1])/m_tbinWid[lay][1] + 0.5);
55 sprintf(hname, "mxt%02d_%02d_%d_mid", lay, iEntr, lr);
56 m_pfMid[lay][iEntr][lr] = new TProfile(hname, hname, xbinM, m_tbinLim[lay][1], m_tbinLim[lay][2]);
57 m_fdPf->Add(m_pfMid[lay][iEntr][lr]);
58
59 int xbinF = (int)((m_tbinLim[lay][3] - m_tbinLim[lay][2])/m_tbinWid[lay][2] + 0.5);
60 sprintf(hname, "mxt%02d_%02d_%d_far", lay, iEntr, lr);
61 m_pfFar[lay][iEntr][lr] = new TProfile(hname, hname, xbinF, m_tbinLim[lay][2], m_tbinLim[lay][3]);
62 m_fdPf->Add(m_pfFar[lay][iEntr][lr]);
63 }
64 }
65 }
66}
67
68void XtInteCalib::mergeHist(TFile* fhist){
69 CalibBase::mergeHist(fhist);
70
71 double tdr, doca;
72 char hname[200];
73 TProfile* pr;
74 TFolder* fd = (TFolder*)fhist->Get("fdProfile");
75 for(int lay=0; lay<NLAYER; lay++){
76 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
77 for(int lr=0; lr<2; lr++){
78 if((m_grXt[lay][iEntr][lr]->GetN()) < m_nMaxGrPoint){
79 sprintf(hname, "xt%02d_%02d_%d_gr", lay, iEntr, lr);
80 TGraph* gr = (TGraph*)fd->FindObjectAny(hname);
81 int nPoint = gr->GetN();
82 for(int i=0; i<nPoint; i++){
83 gr->GetPoint(i, tdr, doca);
84 int np = m_grXt[lay][iEntr][lr]->GetN();
85 m_grXt[lay][iEntr][lr]->SetPoint(np, tdr, doca);
86 }
87 }
88 sprintf(hname, "xt%02d_%02d_%d_near", lay, iEntr, lr);
89 pr = (TProfile*)fd->FindObjectAny(hname);
90 m_pfNear[lay][iEntr][lr]->Add(pr);
91
92 sprintf(hname, "xt%02d_%02d_%d_mid", lay, iEntr, lr);
93 pr = (TProfile*)fd->FindObjectAny(hname);
94 m_pfMid[lay][iEntr][lr]->Add(pr);
95
96 sprintf(hname, "xt%02d_%02d_%d_far", lay, iEntr, lr);
97 pr = (TProfile*)fd->FindObjectAny(hname);
98 m_pfFar[lay][iEntr][lr]->Add(pr);
99 }
100 }
101 }
102}
103
104void XtInteCalib::calib(MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList){
105 CalibBase::calib(calconst, newXtList, r2tList);
106 bool fgOldXt = saveOldXt(newXtList);
107 newXtList->Clear();
108
109 TGraph* grXtfit[NLAYER][NENTRXT][2];
110 char hname[200];
111 for(int lay=0; lay<NLAYER; lay++){
112 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
113 for(int lr=0; lr<2; lr++){
114 m_vt.clear();
115 m_vd.clear();
116 m_entries.clear();
117 for(int iPf=0; iPf<3; iPf++){
118 TProfile* pro;
119 if(0==iPf) pro = m_pfNear[lay][iEntr][lr];
120 else if(1==iPf) pro = m_pfMid[lay][iEntr][lr];
121 else pro = m_pfFar[lay][iEntr][lr];
122
123 int nbin = pro->GetNbinsX();
124 for(int i=0; i<nbin; i++){
125 double tt = pro->GetBinCenter(i+1);
126 double dd = pro->GetBinContent(i+1);
127 double entries = pro->GetBinEntries(i+1);
128 if(entries > 10){
129 m_vt.push_back(tt);
130 m_vd.push_back(dd);
131 m_entries.push_back(entries);
132 }
133 }
134 }
135 unsigned vsize = m_vt.size();
136 if(vsize > 10){
137 for(int i=0; i<2; i++){
138 double slope = (m_vd[vsize-1]-m_vd[vsize-2])/(m_vt[vsize-1]-m_vt[vsize-2]);
139 if(fabs(slope)>0.04){ // 0.8mm/20ns
140 m_vt.pop_back();
141 m_vd.pop_back();
142 m_entries.pop_back();
143 vsize = m_vt.size();
144 }
145 }
146 }
147 sprintf(hname, "grXtFit%02d_%02d_%d", lay, iEntr, lr);
148 grXtfit[lay][iEntr][lr] = new TGraph();
149 grXtfit[lay][iEntr][lr]->SetName(hname);
150 grXtfit[lay][iEntr][lr]->SetMarkerStyle(20);
151 m_fgFit[lay][iEntr][lr] = getXt(lay, iEntr, lr, grXtfit[lay][iEntr][lr]);
152 }
153 }
154 }
155
156 double tdr, doca;
157 for(int lay=0; lay<NLAYER; lay++){
158 double tCut = 500.0;
159 if(lay<8) tCut = 400.0;
160 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
161 for(int lr=0; lr<2; lr++){
162 if(!m_fgFit[lay][iEntr][lr]){
163 int iEntrNew = findXtEntr(lay, iEntr, lr);
164 if(-1 != iEntrNew){
165 int npoint = grXtfit[lay][iEntrNew][lr]->GetN();
166 for(int i=0; i<npoint; i++){
167 grXtfit[lay][iEntrNew][lr]->GetPoint(i, tdr, doca);
168 grXtfit[lay][iEntr][lr]->SetPoint(i, tdr, doca);
169 }
170 } else if(fgOldXt) {
171 cout << grXtfit[lay][iEntr][lr]->GetName() << " use old x-t" << endl;
172 int npoint = m_grXtOld[lay][iEntr][lr]->GetN();
173 for(int i=0; i<npoint; i++){
174 m_grXtOld[lay][iEntr][lr]->GetPoint(i, tdr, doca);
175 grXtfit[lay][iEntr][lr]->SetPoint(i, tdr, doca);
176 }
177 }
178 }
179 int nn = grXtfit[lay][iEntr][lr]->GetN();
180 double tmax, dmax;
181 grXtfit[lay][iEntr][lr]->GetPoint(nn-1, tmax, dmax);
182 if(tmax > tCut) m_fgEdge[lay][iEntr][lr] = true;
183 else m_fgEdge[lay][iEntr][lr] = false;
184 }
185 }
186 }
187
188 for(int lay=0; lay<NLAYER; lay++){
189 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
190 for(int lr=0; lr<2; lr++){
191 if(!m_fgEdge[lay][iEntr][lr]){
192 int iEntrNew = findXtEntrEdge(lay, iEntr, lr);
193 if(-1 != iEntrNew){
194 double t1, d1;
195 int n1 = grXtfit[lay][iEntr][lr]->GetN();
196 grXtfit[lay][iEntr][lr]->GetPoint(n1-1, t1, d1);
197 double t2, d2;
198 int n2 = grXtfit[lay][iEntrNew][lr]->GetN();
199 for(int i=0; i<n2; i++){
200 grXtfit[lay][iEntrNew][lr]->GetPoint(i, t2, d2);
201 if(t2 > t1){
202 grXtfit[lay][iEntr][lr]->SetPoint(n1, t2, d2);
203 n1++;
204 }
205 }
206 }
207 }
208 }
209 }
210 }
211
212 TTree* xttr[NLAYER][NENTRXT][2];
213 for(int lay=0; lay<NLAYER; lay++){
214 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
215 for(int lr=0; lr<2; lr++){
216 sprintf(hname, "trNewXt%02d_%02d_%d", lay, iEntr, lr);
217 xttr[lay][iEntr][lr] = new TTree(hname, hname);
218 xttr[lay][iEntr][lr]->Branch("t", &tdr, "t/D");
219 xttr[lay][iEntr][lr]->Branch("d", &doca, "d/D");
220 if(0 == gFgCalib[lay]){
221 int npoint = m_grXtOld[lay][iEntr][lr]->GetN();
222 for(int i=0; i<npoint; i++){
223 m_grXtOld[lay][iEntr][lr]->GetPoint(i, tdr, doca);
224 xttr[lay][iEntr][lr]->Fill();
225 }
226 } else{
227 int npoint = grXtfit[lay][iEntr][lr]->GetN();
228 for(int i=0; i<npoint; i++){
229 grXtfit[lay][iEntr][lr]->GetPoint(i, tdr, doca);
230 xttr[lay][iEntr][lr]->Fill();
231 }
232 }
233 newXtList->Add(xttr[lay][iEntr][lr]);
234 }
235 }
236 }
237
238 for(int lay=0; lay<NLAYER; lay++){
239 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
240 for(int lr=0; lr<2; lr++){
241 delete grXtfit[lay][iEntr][lr];
242 }
243 }
244 }
245
246// fnewXt.Close();
247 renameHist();
248}
249
250void XtInteCalib::renameHist(){
251 char hname[200];
252 m_fdPf->SetName("fdProfile");
253 for(int lay=0; lay<NLAYER; lay++){
254 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
255 for(int lr=0; lr<2; lr++){
256 sprintf(hname, "xt%02d_%02d_%d_gr", lay, iEntr, lr);
257 m_grXt[lay][iEntr][lr] -> SetName(hname);
258 sprintf(hname, "xt%02d_%02d_%d_near", lay, iEntr, lr);
259 m_pfNear[lay][iEntr][lr] -> SetName(hname);
260 sprintf(hname, "xt%02d_%02d_%d_mid", lay, iEntr, lr);
261 m_pfMid[lay][iEntr][lr] -> SetName(hname);
262 sprintf(hname, "xt%02d_%02d_%d_far", lay, iEntr, lr);
263 m_pfFar[lay][iEntr][lr] -> SetName(hname);
264 }
265 }
266 }
267}
268
269bool XtInteCalib::saveOldXt(TObjArray* newXtList){
270 char hname[200];
271 Int_t entries = newXtList->GetEntries();
272 cout << "entries of newXtList " << entries << endl;
273 if(entries < 1548){
274 cout << "can not get old x-t" << endl;
275 return false;
276 }
277 for(int lay=0; lay<NLAYER; lay++){
278 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
279 for(int lr=0; lr<2; lr++){
280 sprintf(hname, "grXtOld%02d_%02d_%d", lay, iEntr, lr);
281 m_grXtOld[lay][iEntr][lr] = new TGraph();
282 m_grXtOld[lay][iEntr][lr]->SetName(hname);
283
284 double t, d;
285 sprintf(hname, "trNewXt%02d_%02d_%d", lay, iEntr, lr);
286 TTree* tr = (TTree*)(newXtList->FindObject(hname));
287// cout << setw(15) << tr->GetName() << setw(15) << tr->GetEntries() << endl;
288 tr -> SetBranchAddress("t", &t);
289 tr -> SetBranchAddress("d", &d);
290 int nPoint = tr -> GetEntries();
291 if(nPoint < 20){
292 cout << "can not get old x-t: " << hname << endl;
293 return false;
294 }
295 for(int i=0; i<nPoint; i++){
296 tr->GetEntry(i);
297 m_grXtOld[lay][iEntr][lr]->SetPoint(i, t, d);
298 }
299 }
300 }
301 }
302// for(int lay=0; lay<NLAYER; lay++){
303// for(int iEntr=0; iEntr<NENTRXT; iEntr++){
304// for(int lr=0; lr<2; lr++){
305// sprintf(hname, "trNewXt%02d_%02d_%d", lay, iEntr, lr);
306// delete (newXtList->FindObject(hname));
307// }
308// }
309// }
310 return true;
311}
312
313bool XtInteCalib::getXt(int lay, int iEntr, int lr, TGraph* gr){
314 unsigned vsize = m_vt.size();
315 if(vsize < 15) return false;
316
317 double tmax = m_vt[vsize-1];
318 double tm0 = 300.0; // 280->300
319 double tm1 = 500.0;
320 double tm2 = 700.0;
321 if(lay<8){
322 tm0 = 200.0; // 180->200
323 tm1 = 300.0;
324 tm2 = 540.0;
325 }
326 int n0 = 0;
327 for(unsigned i=0; i<vsize; i++){
328 if(m_vt[i] < tm0) n0++;
329 }
330
331 int nCut = 30; // 25
332 if(lay<8) nCut = 20; // 15
333
334 double entries1 = 0.0;
335 double entries2 = 0.0;
336 for(unsigned i=0; i<vsize; i++){
337 entries1 += m_entries[i];
338 if(m_vt[i] < 150.) entries2 += m_entries[i];
339 }
340 if((entries1*0.9) < entries2) return false;
341
342 if(n0 < nCut) return false;
343 if(tmax < tm1) return funXt0(lay,iEntr,lr, gr);
344 else return funXt1(lay,iEntr,lr, gr);
345}
346
347bool XtInteCalib::funXt0(int lay, int iEntr, int lr, TGraph* gr){
348 double tCut = 300.0;
349 if(lay<8) tCut = 200.0;
350
351 int npoint = 0;
352 if(m_vt[0] > 0.0){
353 gr->SetPoint(0, 0, 0);
354 npoint++;
355 }
356 for(unsigned i=0; i<m_vt.size(); i++){
357 double delt = 0.0;
358 if(i>10) delt = m_vt[i] - m_vt[i-1];
359 if((i>10) && ((delt>100.) || ((delt>60.) && (m_vt[i-1]<tCut)))) break;
360
361 if(m_vt[i] < 0.0) gr->SetPoint(npoint, m_vt[i], 0.0);
362 else gr->SetPoint(npoint, m_vt[i], m_vd[i]);
363 npoint++;
364 }
365 return true;
366}
367
368bool XtInteCalib::funXt1(int lay, int iEntr, int lr, TGraph* gr){
369 unsigned vsize = m_vt.size();
370 double tmax = m_vt[vsize-1];
371 double t1;
372 if(lay<8){
373 if(tmax<540) t1 = 300.;
374 else t1 = 540.;
375 }else{
376 if(tmax<700) t1 = 500.;
377 else t1 = 660.;
378 }
379
380 bool fgfit = false;
381 int np = 0;
382 double c0 = 0.0;
383 double c1 = 0.0;
384 TGraph* grPol1 = new TGraph();
385 for(unsigned i=0; i<vsize; i++){
386 if(m_vt[i] >= t1){
387 grPol1->SetPoint(np, m_vt[i], m_vd[i]);
388 np++;
389 }
390 }
391 double t2;
392 if(np<5){
393 t2 = t1;
394 } else{
395 double x2;
396 grPol1->GetPoint(0,t2,x2);
397 grPol1->Fit("pol1","Q","",t2, m_vt[vsize-1]);
398 c0 = grPol1->GetFunction("pol1")->GetParameter(0);
399 c1 = grPol1->GetFunction("pol1")->GetParameter(1);
400
401 if(c1<0){
402 grPol1->Fit("pol0","Q","",t2, m_vt[vsize-1]);
403 c0 = grPol1->GetFunction("pol0")->GetParameter(0);
404 c1 = 0.0;
405 }
406 fgfit = true;
407 }
408
409 double tCut = 300.0;
410 if(lay<8) tCut = 200.0;
411
412 int npoint = 0;
413 if(m_vt[0] > 0.0){
414 gr->SetPoint(0, 0, 0);
415 npoint++;
416 }
417 for(unsigned i=0; i<vsize; i++){
418 double delt = 0.0;
419 if(i>10) delt = m_vt[i] - m_vt[i-1];
420 if((i>10) && ((delt>100.) || ((delt>60.) && (m_vt[i-1]<tCut)))) break;
421
422 if(m_vt[i] < t2){
423 if(m_vt[i] < 0.0) gr->SetPoint(npoint, m_vt[i], 0.0);
424 else gr->SetPoint(npoint, m_vt[i], m_vd[i]);
425 npoint++;
426 } else{
427 double dist;
428 if(!fgfit) dist = m_vd[i];
429 else dist = c1*m_vt[i] + c0;
430 gr->SetPoint(npoint, m_vt[i], dist);
431 npoint++;
432 }
433 }
434 return true;
435}
436
437
438int XtInteCalib::findXtEntr(int lay, int iEntr, int lr) const {
439 int id0 = 8;
440 int id1 = 9;
441 int idmax = 17;
442 int entrId = -1;
443 if(iEntr <= id0){
444 int id = -1;
445 for(int i=iEntr; i<=id0; i++){
446 if(m_fgFit[lay][i][lr]){
447 id = i;
448 break;
449 }
450 }
451 if(-1 != id) entrId = id;
452 else{
453 for(int i=iEntr; i>=0; i--){
454 if(m_fgFit[lay][i][lr]){
455 id = i;
456 break;
457 }
458 }
459 if(-1 != id) entrId = id;
460 else{
461 for(int i=id1; i<=idmax; i++){
462 if(m_fgFit[lay][i][lr]){
463 id = i;
464 break;
465 }
466 }
467 entrId = id;
468 }
469 }
470 } else{
471 int id = -1;
472 for(int i=iEntr; i>=id1; i--){
473 if(m_fgFit[lay][i][lr]){
474 id = i;
475 break;
476 }
477 }
478 if(-1 != id) entrId = id;
479 else{
480 for(int i=iEntr; i<idmax; i++){
481 if(m_fgFit[lay][i][lr]){
482 id = i;
483 break;
484 }
485 }
486 if(-1 != id) entrId = id;
487 else{
488 for(int i=id1; i>=0; i--){
489 if(m_fgFit[lay][i][lr]){
490 id = i;
491 break;
492 }
493 }
494 entrId = id;
495 }
496 }
497 }
498 if(-1 == entrId){
499 cout << "find EntrId error " << "layer " << lay << " iEntr " << iEntr << " lr " << lr << endl;
500 }
501
502 return entrId;
503}
504
505int XtInteCalib::findXtEntrEdge(int lay, int iEntr, int lr) const {
506 int id0 = 8;
507 int id1 = 9;
508 int idmax = 17;
509 int entrId = -1;
510 if(iEntr <= id0){
511 int id = -1;
512 for(int i=iEntr; i<=id0; i++){
513 if(m_fgEdge[lay][i][lr]){
514 id = i;
515 break;
516 }
517 }
518 if(-1 != id) entrId = id;
519 else{
520 for(int i=iEntr; i>=0; i--){
521 if(m_fgEdge[lay][i][lr]){
522 id = i;
523 break;
524 }
525 }
526 if(-1 != id) entrId = id;
527 else{
528 for(int i=id1; i<=idmax; i++){
529 if(m_fgEdge[lay][i][lr]){
530 id = i;
531 break;
532 }
533 }
534 entrId = id;
535 }
536 }
537 } else{
538 int id = -1;
539 for(int i=iEntr; i>=id1; i--){
540 if(m_fgEdge[lay][i][lr]){
541 id = i;
542 break;
543 }
544 }
545 if(-1 != id) entrId = id;
546 else{
547 for(int i=iEntr; i<idmax; i++){
548 if(m_fgEdge[lay][i][lr]){
549 id = i;
550 break;
551 }
552 }
553 if(-1 != id) entrId = id;
554 else{
555 for(int i=id1; i>=0; i--){
556 if(m_fgEdge[lay][i][lr]){
557 id = i;
558 break;
559 }
560 }
561 entrId = id;
562 }
563 }
564 }
565 if(-1 == entrId){
566 cout << "find EntrId error for cell edge " << "layer " << lay << " iEntr " << iEntr << " lr " << lr << endl;
567 }
568
569 return entrId;
570}
c1
TCanvas * c1
Definition CgemLUTWriter.cxx:75
SetBranchAddress
data SetBranchAddress("time",&time)
gr
TGraph * gr
Definition DataBase/tau_mode.c:69
n2
int n2
Definition SD0Tag.cxx:55
n1
int n1
Definition SD0Tag.cxx:54
CalibBase::init
virtual void init(TObjArray *hlist, MdcCosGeom *pGeom)=0
Definition CalibBase.cpp:14
CalibBase::mergeHist
virtual void mergeHist(TFile *fhist)=0
Definition CalibBase.cpp:365
CalibBase::calib
virtual void calib(MdcCalibConst *calconst, TObjArray *newXtList, TObjArray *r2tList)=0
Definition CalibBase.cpp:711
XtInteCalib::calib
void calib(MdcCalibConst *calconst, TObjArray *newXtList, TObjArray *r2tList)
Definition XtInteCalib.cpp:104
XtInteCalib::mergeHist
void mergeHist(TFile *fhist)
Definition XtInteCalib.cpp:68
XtInteCalib::init
void init(TObjArray *hlist, MdcCosGeom *pGeom)
Definition XtInteCalib.cpp:29
t
int t()
Definition t.c:1