00001 #include "include/XtInteCalib.h"
00002 #include "include/fun.h"
00003
00004 #include <cmath>
00005
00006 #include "TMinuit.h"
00007 #include "TF1.h"
00008 #include "TTree.h"
00009
00010 XtInteCalib::XtInteCalib(){
00011 m_nMaxGrPoint = 5000;
00012 for(int lay=0; lay<NLAYER; lay++){
00013 m_tbinWid[lay][0] = 5.0;
00014 m_tbinWid[lay][1] = 10.0;
00015 m_tbinWid[lay][2] = 20.0;
00016
00017 m_tbinLim[lay][0] = -10.0;
00018 m_tbinLim[lay][1] = 30.0;
00019 if(lay < 8) m_tbinLim[lay][2] = 210.0;
00020 else m_tbinLim[lay][2] = 350.0;
00021 m_tbinLim[lay][3] = 990.0;
00022 }
00023 cout << "Calibration type: XtInteCalib" << endl;
00024 }
00025
00026 XtInteCalib::~XtInteCalib(){
00027 }
00028
00029 void XtInteCalib::init(TObjArray* hlist, MdcCosGeom* pGeom){
00030 CalibBase::init(hlist, pGeom);
00031
00032 m_fdPf = new TFolder("mfdProfile", "fdProfile");
00033 hlist -> Add(m_fdPf);
00034
00035 m_haxis = new TH2F("axis", "", 200, -50, 1000, 50, 0, 10);
00036 m_haxis -> SetStats(0);
00037 m_fdPf -> Add(m_haxis);
00038
00039 char hname[200];
00040 for(int lay=0; lay<NLAYER; lay++){
00041 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00042 for(int lr=0; lr<2; lr++){
00043 sprintf(hname, "mxt%02d_%02d_%d_gr", lay, iEntr, lr);
00044 m_grXt[lay][iEntr][lr] = new TGraph();
00045 m_grXt[lay][iEntr][lr]->SetName(hname);
00046 m_grXt[lay][iEntr][lr]->SetMarkerColor(2);
00047 m_fdPf->Add(m_grXt[lay][iEntr][lr]);
00048
00049 int xbinN = (int)((m_tbinLim[lay][1] - m_tbinLim[lay][0])/m_tbinWid[lay][0] + 0.5);
00050 sprintf(hname, "mxt%02d_%02d_%d_near", lay, iEntr, lr);
00051 m_pfNear[lay][iEntr][lr] = new TProfile(hname, hname, xbinN, m_tbinLim[lay][0], m_tbinLim[lay][1]);
00052 m_fdPf->Add(m_pfNear[lay][iEntr][lr]);
00053
00054 int xbinM = (int)((m_tbinLim[lay][2] - m_tbinLim[lay][1])/m_tbinWid[lay][1] + 0.5);
00055 sprintf(hname, "mxt%02d_%02d_%d_mid", lay, iEntr, lr);
00056 m_pfMid[lay][iEntr][lr] = new TProfile(hname, hname, xbinM, m_tbinLim[lay][1], m_tbinLim[lay][2]);
00057 m_fdPf->Add(m_pfMid[lay][iEntr][lr]);
00058
00059 int xbinF = (int)((m_tbinLim[lay][3] - m_tbinLim[lay][2])/m_tbinWid[lay][2] + 0.5);
00060 sprintf(hname, "mxt%02d_%02d_%d_far", lay, iEntr, lr);
00061 m_pfFar[lay][iEntr][lr] = new TProfile(hname, hname, xbinF, m_tbinLim[lay][2], m_tbinLim[lay][3]);
00062 m_fdPf->Add(m_pfFar[lay][iEntr][lr]);
00063 }
00064 }
00065 }
00066 }
00067
00068 void XtInteCalib::mergeHist(TFile* fhist){
00069 CalibBase::mergeHist(fhist);
00070
00071 double tdr, doca;
00072 char hname[200];
00073 TProfile* pr;
00074 TFolder* fd = (TFolder*)fhist->Get("fdProfile");
00075 for(int lay=0; lay<NLAYER; lay++){
00076 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00077 for(int lr=0; lr<2; lr++){
00078 if((m_grXt[lay][iEntr][lr]->GetN()) < m_nMaxGrPoint){
00079 sprintf(hname, "xt%02d_%02d_%d_gr", lay, iEntr, lr);
00080 TGraph* gr = (TGraph*)fd->FindObjectAny(hname);
00081 int nPoint = gr->GetN();
00082 for(int i=0; i<nPoint; i++){
00083 gr->GetPoint(i, tdr, doca);
00084 int np = m_grXt[lay][iEntr][lr]->GetN();
00085 m_grXt[lay][iEntr][lr]->SetPoint(np, tdr, doca);
00086 }
00087 }
00088 sprintf(hname, "xt%02d_%02d_%d_near", lay, iEntr, lr);
00089 pr = (TProfile*)fd->FindObjectAny(hname);
00090 m_pfNear[lay][iEntr][lr]->Add(pr);
00091
00092 sprintf(hname, "xt%02d_%02d_%d_mid", lay, iEntr, lr);
00093 pr = (TProfile*)fd->FindObjectAny(hname);
00094 m_pfMid[lay][iEntr][lr]->Add(pr);
00095
00096 sprintf(hname, "xt%02d_%02d_%d_far", lay, iEntr, lr);
00097 pr = (TProfile*)fd->FindObjectAny(hname);
00098 m_pfFar[lay][iEntr][lr]->Add(pr);
00099 }
00100 }
00101 }
00102 }
00103
00104 void XtInteCalib::calib(MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList){
00105 CalibBase::calib(calconst, newXtList, r2tList);
00106 bool fgOldXt = saveOldXt(newXtList);
00107 newXtList->Clear();
00108
00109 TGraph* grXtfit[NLAYER][NENTRXT][2];
00110 char hname[200];
00111 for(int lay=0; lay<NLAYER; lay++){
00112 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00113 for(int lr=0; lr<2; lr++){
00114 m_vt.clear();
00115 m_vd.clear();
00116 m_entries.clear();
00117 for(int iPf=0; iPf<3; iPf++){
00118 TProfile* pro;
00119 if(0==iPf) pro = m_pfNear[lay][iEntr][lr];
00120 else if(1==iPf) pro = m_pfMid[lay][iEntr][lr];
00121 else pro = m_pfFar[lay][iEntr][lr];
00122
00123 int nbin = pro->GetNbinsX();
00124 for(int i=0; i<nbin; i++){
00125 double tt = pro->GetBinCenter(i+1);
00126 double dd = pro->GetBinContent(i+1);
00127 double entries = pro->GetBinEntries(i+1);
00128 if(entries > 10){
00129 m_vt.push_back(tt);
00130 m_vd.push_back(dd);
00131 m_entries.push_back(entries);
00132 }
00133 }
00134 }
00135 unsigned vsize = m_vt.size();
00136 if(vsize > 10){
00137 for(int i=0; i<2; i++){
00138 double slope = (m_vd[vsize-1]-m_vd[vsize-2])/(m_vt[vsize-1]-m_vt[vsize-2]);
00139 if(fabs(slope)>0.04){
00140 m_vt.pop_back();
00141 m_vd.pop_back();
00142 m_entries.pop_back();
00143 vsize = m_vt.size();
00144 }
00145 }
00146 }
00147 sprintf(hname, "grXtFit%02d_%02d_%d", lay, iEntr, lr);
00148 grXtfit[lay][iEntr][lr] = new TGraph();
00149 grXtfit[lay][iEntr][lr]->SetName(hname);
00150 grXtfit[lay][iEntr][lr]->SetMarkerStyle(20);
00151 m_fgFit[lay][iEntr][lr] = getXt(lay, iEntr, lr, grXtfit[lay][iEntr][lr]);
00152 }
00153 }
00154 }
00155
00156 double tdr, doca;
00157 for(int lay=0; lay<NLAYER; lay++){
00158 double tCut = 500.0;
00159 if(lay<8) tCut = 400.0;
00160 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00161 for(int lr=0; lr<2; lr++){
00162 if(!m_fgFit[lay][iEntr][lr]){
00163 int iEntrNew = findXtEntr(lay, iEntr, lr);
00164 if(-1 != iEntrNew){
00165 int npoint = grXtfit[lay][iEntrNew][lr]->GetN();
00166 for(int i=0; i<npoint; i++){
00167 grXtfit[lay][iEntrNew][lr]->GetPoint(i, tdr, doca);
00168 grXtfit[lay][iEntr][lr]->SetPoint(i, tdr, doca);
00169 }
00170 } else if(fgOldXt) {
00171 cout << grXtfit[lay][iEntr][lr]->GetName() << " use old x-t" << endl;
00172 int npoint = m_grXtOld[lay][iEntr][lr]->GetN();
00173 for(int i=0; i<npoint; i++){
00174 m_grXtOld[lay][iEntr][lr]->GetPoint(i, tdr, doca);
00175 grXtfit[lay][iEntr][lr]->SetPoint(i, tdr, doca);
00176 }
00177 }
00178 }
00179 int nn = grXtfit[lay][iEntr][lr]->GetN();
00180 double tmax, dmax;
00181 grXtfit[lay][iEntr][lr]->GetPoint(nn-1, tmax, dmax);
00182 if(tmax > tCut) m_fgEdge[lay][iEntr][lr] = true;
00183 else m_fgEdge[lay][iEntr][lr] = false;
00184 }
00185 }
00186 }
00187
00188 for(int lay=0; lay<NLAYER; lay++){
00189 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00190 for(int lr=0; lr<2; lr++){
00191 if(!m_fgEdge[lay][iEntr][lr]){
00192 int iEntrNew = findXtEntrEdge(lay, iEntr, lr);
00193 if(-1 != iEntrNew){
00194 double t1, d1;
00195 int n1 = grXtfit[lay][iEntr][lr]->GetN();
00196 grXtfit[lay][iEntr][lr]->GetPoint(n1-1, t1, d1);
00197 double t2, d2;
00198 int n2 = grXtfit[lay][iEntrNew][lr]->GetN();
00199 for(int i=0; i<n2; i++){
00200 grXtfit[lay][iEntrNew][lr]->GetPoint(i, t2, d2);
00201 if(t2 > t1){
00202 grXtfit[lay][iEntr][lr]->SetPoint(n1, t2, d2);
00203 n1++;
00204 }
00205 }
00206 }
00207 }
00208 }
00209 }
00210 }
00211
00212 TTree* xttr[NLAYER][NENTRXT][2];
00213 for(int lay=0; lay<NLAYER; lay++){
00214 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00215 for(int lr=0; lr<2; lr++){
00216 sprintf(hname, "trNewXt%02d_%02d_%d", lay, iEntr, lr);
00217 xttr[lay][iEntr][lr] = new TTree(hname, hname);
00218 xttr[lay][iEntr][lr]->Branch("t", &tdr, "t/D");
00219 xttr[lay][iEntr][lr]->Branch("d", &doca, "d/D");
00220 if(0 == gFgCalib[lay]){
00221 int npoint = m_grXtOld[lay][iEntr][lr]->GetN();
00222 for(int i=0; i<npoint; i++){
00223 m_grXtOld[lay][iEntr][lr]->GetPoint(i, tdr, doca);
00224 xttr[lay][iEntr][lr]->Fill();
00225 }
00226 } else{
00227 int npoint = grXtfit[lay][iEntr][lr]->GetN();
00228 for(int i=0; i<npoint; i++){
00229 grXtfit[lay][iEntr][lr]->GetPoint(i, tdr, doca);
00230 xttr[lay][iEntr][lr]->Fill();
00231 }
00232 }
00233 newXtList->Add(xttr[lay][iEntr][lr]);
00234 }
00235 }
00236 }
00237
00238 for(int lay=0; lay<NLAYER; lay++){
00239 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00240 for(int lr=0; lr<2; lr++){
00241 delete grXtfit[lay][iEntr][lr];
00242 }
00243 }
00244 }
00245
00246
00247 renameHist();
00248 }
00249
00250 void XtInteCalib::renameHist(){
00251 char hname[200];
00252 m_fdPf->SetName("fdProfile");
00253 for(int lay=0; lay<NLAYER; lay++){
00254 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00255 for(int lr=0; lr<2; lr++){
00256 sprintf(hname, "xt%02d_%02d_%d_gr", lay, iEntr, lr);
00257 m_grXt[lay][iEntr][lr] -> SetName(hname);
00258 sprintf(hname, "xt%02d_%02d_%d_near", lay, iEntr, lr);
00259 m_pfNear[lay][iEntr][lr] -> SetName(hname);
00260 sprintf(hname, "xt%02d_%02d_%d_mid", lay, iEntr, lr);
00261 m_pfMid[lay][iEntr][lr] -> SetName(hname);
00262 sprintf(hname, "xt%02d_%02d_%d_far", lay, iEntr, lr);
00263 m_pfFar[lay][iEntr][lr] -> SetName(hname);
00264 }
00265 }
00266 }
00267 }
00268
00269 bool XtInteCalib::saveOldXt(TObjArray* newXtList){
00270 char hname[200];
00271 Int_t entries = newXtList->GetEntries();
00272 cout << "entries of newXtList " << entries << endl;
00273 if(entries < 1548){
00274 cout << "can not get old x-t" << endl;
00275 return false;
00276 }
00277 for(int lay=0; lay<NLAYER; lay++){
00278 for(int iEntr=0; iEntr<NENTRXT; iEntr++){
00279 for(int lr=0; lr<2; lr++){
00280 sprintf(hname, "grXtOld%02d_%02d_%d", lay, iEntr, lr);
00281 m_grXtOld[lay][iEntr][lr] = new TGraph();
00282 m_grXtOld[lay][iEntr][lr]->SetName(hname);
00283
00284 double t, d;
00285 sprintf(hname, "trNewXt%02d_%02d_%d", lay, iEntr, lr);
00286 TTree* tr = (TTree*)(newXtList->FindObject(hname));
00287
00288 tr -> SetBranchAddress("t", &t);
00289 tr -> SetBranchAddress("d", &d);
00290 int nPoint = tr -> GetEntries();
00291 if(nPoint < 20){
00292 cout << "can not get old x-t: " << hname << endl;
00293 return false;
00294 }
00295 for(int i=0; i<nPoint; i++){
00296 tr->GetEntry(i);
00297 m_grXtOld[lay][iEntr][lr]->SetPoint(i, t, d);
00298 }
00299 }
00300 }
00301 }
00302
00303
00304
00305
00306
00307
00308
00309
00310 return true;
00311 }
00312
00313 bool XtInteCalib::getXt(int lay, int iEntr, int lr, TGraph* gr){
00314 unsigned vsize = m_vt.size();
00315 if(vsize < 15) return false;
00316
00317 double tmax = m_vt[vsize-1];
00318 double tm0 = 300.0;
00319 double tm1 = 500.0;
00320 double tm2 = 700.0;
00321 if(lay<8){
00322 tm0 = 200.0;
00323 tm1 = 300.0;
00324 tm2 = 540.0;
00325 }
00326 int n0 = 0;
00327 for(unsigned i=0; i<vsize; i++){
00328 if(m_vt[i] < tm0) n0++;
00329 }
00330
00331 int nCut = 30;
00332 if(lay<8) nCut = 20;
00333
00334 double entries1 = 0.0;
00335 double entries2 = 0.0;
00336 for(unsigned i=0; i<vsize; i++){
00337 entries1 += m_entries[i];
00338 if(m_vt[i] < 150.) entries2 += m_entries[i];
00339 }
00340 if((entries1*0.9) < entries2) return false;
00341
00342 if(n0 < nCut) return false;
00343 if(tmax < tm1) return funXt0(lay,iEntr,lr, gr);
00344 else return funXt1(lay,iEntr,lr, gr);
00345 }
00346
00347 bool XtInteCalib::funXt0(int lay, int iEntr, int lr, TGraph* gr){
00348 double tCut = 300.0;
00349 if(lay<8) tCut = 200.0;
00350
00351 int npoint = 0;
00352 if(m_vt[0] > 0.0){
00353 gr->SetPoint(0, 0, 0);
00354 npoint++;
00355 }
00356 for(unsigned i=0; i<m_vt.size(); i++){
00357 double delt = 0.0;
00358 if(i>10) delt = m_vt[i] - m_vt[i-1];
00359 if((i>10) && ((delt>100.) || ((delt>60.) && (m_vt[i-1]<tCut)))) break;
00360
00361 if(m_vt[i] < 0.0) gr->SetPoint(npoint, m_vt[i], 0.0);
00362 else gr->SetPoint(npoint, m_vt[i], m_vd[i]);
00363 npoint++;
00364 }
00365 return true;
00366 }
00367
00368 bool XtInteCalib::funXt1(int lay, int iEntr, int lr, TGraph* gr){
00369 unsigned vsize = m_vt.size();
00370 double tmax = m_vt[vsize-1];
00371 double t1;
00372 if(lay<8){
00373 if(tmax<540) t1 = 300.;
00374 else t1 = 540.;
00375 }else{
00376 if(tmax<700) t1 = 500.;
00377 else t1 = 660.;
00378 }
00379
00380 bool fgfit = false;
00381 int np = 0;
00382 double c0 = 0.0;
00383 double c1 = 0.0;
00384 TGraph* grPol1 = new TGraph();
00385 for(unsigned i=0; i<vsize; i++){
00386 if(m_vt[i] >= t1){
00387 grPol1->SetPoint(np, m_vt[i], m_vd[i]);
00388 np++;
00389 }
00390 }
00391 double t2;
00392 if(np<5){
00393 t2 = t1;
00394 } else{
00395 double x2;
00396 grPol1->GetPoint(0,t2,x2);
00397 grPol1->Fit("pol1","Q","",t2, m_vt[vsize-1]);
00398 c0 = grPol1->GetFunction("pol1")->GetParameter(0);
00399 c1 = grPol1->GetFunction("pol1")->GetParameter(1);
00400
00401 if(c1<0){
00402 grPol1->Fit("pol0","Q","",t2, m_vt[vsize-1]);
00403 c0 = grPol1->GetFunction("pol0")->GetParameter(0);
00404 c1 = 0.0;
00405 }
00406 fgfit = true;
00407 }
00408
00409 double tCut = 300.0;
00410 if(lay<8) tCut = 200.0;
00411
00412 int npoint = 0;
00413 if(m_vt[0] > 0.0){
00414 gr->SetPoint(0, 0, 0);
00415 npoint++;
00416 }
00417 for(unsigned i=0; i<vsize; i++){
00418 double delt = 0.0;
00419 if(i>10) delt = m_vt[i] - m_vt[i-1];
00420 if((i>10) && ((delt>100.) || ((delt>60.) && (m_vt[i-1]<tCut)))) break;
00421
00422 if(m_vt[i] < t2){
00423 if(m_vt[i] < 0.0) gr->SetPoint(npoint, m_vt[i], 0.0);
00424 else gr->SetPoint(npoint, m_vt[i], m_vd[i]);
00425 npoint++;
00426 } else{
00427 double dist;
00428 if(!fgfit) dist = m_vd[i];
00429 else dist = c1*m_vt[i] + c0;
00430 gr->SetPoint(npoint, m_vt[i], dist);
00431 npoint++;
00432 }
00433 }
00434 return true;
00435 }
00436
00437
00438 int XtInteCalib::findXtEntr(int lay, int iEntr, int lr) const {
00439 int id0 = 8;
00440 int id1 = 9;
00441 int idmax = 17;
00442 int entrId = -1;
00443 if(iEntr <= id0){
00444 int id = -1;
00445 for(int i=iEntr; i<=id0; i++){
00446 if(m_fgFit[lay][i][lr]){
00447 id = i;
00448 break;
00449 }
00450 }
00451 if(-1 != id) entrId = id;
00452 else{
00453 for(int i=iEntr; i>=0; i--){
00454 if(m_fgFit[lay][i][lr]){
00455 id = i;
00456 break;
00457 }
00458 }
00459 if(-1 != id) entrId = id;
00460 else{
00461 for(int i=id1; i<=idmax; i++){
00462 if(m_fgFit[lay][i][lr]){
00463 id = i;
00464 break;
00465 }
00466 }
00467 entrId = id;
00468 }
00469 }
00470 } else{
00471 int id = -1;
00472 for(int i=iEntr; i>=id1; i--){
00473 if(m_fgFit[lay][i][lr]){
00474 id = i;
00475 break;
00476 }
00477 }
00478 if(-1 != id) entrId = id;
00479 else{
00480 for(int i=iEntr; i<idmax; i++){
00481 if(m_fgFit[lay][i][lr]){
00482 id = i;
00483 break;
00484 }
00485 }
00486 if(-1 != id) entrId = id;
00487 else{
00488 for(int i=id1; i>=0; i--){
00489 if(m_fgFit[lay][i][lr]){
00490 id = i;
00491 break;
00492 }
00493 }
00494 entrId = id;
00495 }
00496 }
00497 }
00498 if(-1 == entrId){
00499 cout << "find EntrId error " << "layer " << lay << " iEntr " << iEntr << " lr " << lr << endl;
00500 }
00501
00502 return entrId;
00503 }
00504
00505 int XtInteCalib::findXtEntrEdge(int lay, int iEntr, int lr) const {
00506 int id0 = 8;
00507 int id1 = 9;
00508 int idmax = 17;
00509 int entrId = -1;
00510 if(iEntr <= id0){
00511 int id = -1;
00512 for(int i=iEntr; i<=id0; i++){
00513 if(m_fgEdge[lay][i][lr]){
00514 id = i;
00515 break;
00516 }
00517 }
00518 if(-1 != id) entrId = id;
00519 else{
00520 for(int i=iEntr; i>=0; i--){
00521 if(m_fgEdge[lay][i][lr]){
00522 id = i;
00523 break;
00524 }
00525 }
00526 if(-1 != id) entrId = id;
00527 else{
00528 for(int i=id1; i<=idmax; i++){
00529 if(m_fgEdge[lay][i][lr]){
00530 id = i;
00531 break;
00532 }
00533 }
00534 entrId = id;
00535 }
00536 }
00537 } else{
00538 int id = -1;
00539 for(int i=iEntr; i>=id1; i--){
00540 if(m_fgEdge[lay][i][lr]){
00541 id = i;
00542 break;
00543 }
00544 }
00545 if(-1 != id) entrId = id;
00546 else{
00547 for(int i=iEntr; i<idmax; i++){
00548 if(m_fgEdge[lay][i][lr]){
00549 id = i;
00550 break;
00551 }
00552 }
00553 if(-1 != id) entrId = id;
00554 else{
00555 for(int i=id1; i>=0; i--){
00556 if(m_fgEdge[lay][i][lr]){
00557 id = i;
00558 break;
00559 }
00560 }
00561 entrId = id;
00562 }
00563 }
00564 }
00565 if(-1 == entrId){
00566 cout << "find EntrId error for cell edge " << "layer " << lay << " iEntr " << iEntr << " lr " << lr << endl;
00567 }
00568
00569 return entrId;
00570 }