00001 #include "include/XtCalib.h"
00002 #include "include/fun.h"
00003
00004 #include "TMinuit.h"
00005
00006 XtCalib::XtCalib(){
00007 cout << "Calibration type: XtCalib" << endl;
00008 }
00009
00010 XtCalib::~XtCalib(){
00011 }
00012
00013 void XtCalib::init(TObjArray* hlist, MdcCosGeom* pGeom){
00014 CalibBase::init(hlist, pGeom);
00015
00016 m_fdXt = new TFolder("mfdxt","fdxt");
00017 hlist->Add(m_fdXt);
00018
00019 char hname[200];
00020 for(int lay=0; lay<43; lay++){
00021 for(int iEntr=0; iEntr<gNEntr[lay]; iEntr++){
00022 for(int lr=0; lr<NLR; lr++){
00023 for(int bin=0; bin<=NXTBIN; bin++){
00024 sprintf(hname, "mHxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, lr, bin);
00025 m_hxt[lay][iEntr][lr][bin] = new TH1D(hname, "", 300, -1.5, 1.5);
00026 m_fdXt->Add(m_hxt[lay][iEntr][lr][bin]);
00027 }
00028 }
00029 }
00030 }
00031 }
00032
00033 void XtCalib::mergeHist(TFile* fhist){
00034 CalibBase::mergeHist(fhist);
00035
00036 char hname[200];
00037 TFolder* fd = (TFolder*)fhist->Get("fdXt");
00038 for(int lay=0; lay<43; lay++){
00039 for(int iEntr=0; iEntr<gNEntr[lay]; iEntr++){
00040 for(int lr=0; lr<NLR; lr++){
00041 for(int bin=0; bin<=NXTBIN; bin++){
00042 sprintf(hname, "Hxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, lr, bin);
00043 TH1F* hist = (TH1F*)fd->FindObjectAny(hname);
00044 m_hxt[lay][iEntr][lr][bin]->Add(hist);
00045
00046
00047
00048
00049 }
00050 }
00051 }
00052 }
00053 }
00054
00055 void XtCalib::calib(MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList){
00056 CalibBase::calib(calconst, newXtList, r2tList);
00057
00058 Int_t ierflg;
00059 Int_t istat;
00060 Int_t nvpar;
00061 Int_t nparx;
00062 Double_t fmin;
00063 Double_t edm;
00064 Double_t errdef;
00065 Double_t arglist[10];
00066
00067 TMinuit* gmxt = new TMinuit(6);
00068 gmxt -> SetPrintLevel(-1);
00069 gmxt -> SetFCN(fcnXT);
00070 gmxt -> SetErrorDef(1.0);
00071 gmxt -> mnparm(0, "xtpar0", 0, 0.1, 0, 0, ierflg);
00072 gmxt -> mnparm(1, "xtpar1", 0, 0.1, 0, 0, ierflg);
00073 gmxt -> mnparm(2, "xtpar2", 0, 0.1, 0, 0, ierflg);
00074 gmxt -> mnparm(3, "xtpar3", 0, 0.1, 0, 0, ierflg);
00075 gmxt -> mnparm(4, "xtpar4", 0, 0.1, 0, 0, ierflg);
00076 gmxt -> mnparm(5, "xtpar5", 0, 0.1, 0, 0, ierflg);
00077 arglist[0] = 0;
00078 gmxt -> mnexcm("SET NOW", arglist, 0, ierflg);
00079
00080 TMinuit* gmxtEd = new TMinuit(1);
00081 gmxtEd -> SetPrintLevel(-1);
00082 gmxtEd -> SetFCN(fcnXtEdge);
00083 gmxtEd -> SetErrorDef(1.0);
00084 gmxtEd -> mnparm(0, "xtpar0", 0, 0.1, 0, 0, ierflg);
00085 arglist[0] = 0;
00086 gmxtEd -> mnexcm("SET NOW", arglist, 0, ierflg);
00087
00088
00089 int i;
00090 Stat_t histEntry;
00091 double xtpar;
00092 double xterr;
00093 double tbcen;
00094 double deltx;
00095 double xcor;
00096 double xerr;
00097 double xtini[8];
00098 double xtfit[8];
00099 ofstream fxtlog("xtlog");
00100 for(int lay=0; lay<43; lay++){
00101 if(0 == gFgCalib[lay]) continue;
00102 for(int iEntr=0; iEntr<gNEntr[lay]; iEntr++){
00103 for(int iLR=0; iLR<NLR; iLR++){
00104 fxtlog << "Layer " << setw(3) << lay << setw(3) << iEntr
00105 << setw(3) << iLR << endl;
00106 for(int ord=0; ord<NXTPAR; ord++){
00107 if(0 == iEntr) xtpar = calconst -> getXtpar(lay, 8, iLR, ord);
00108 else if(1 == iEntr) xtpar = calconst -> getXtpar(lay, 9, iLR, ord);
00109 xtini[ord] = xtpar;
00110 xtfit[ord] = xtpar;
00111 }
00112 Tmax = xtini[6];
00113
00114 for(int bin=0; bin<=NXTBIN; bin++){
00115 histEntry = (int)(m_hxt[lay][iEntr][iLR][bin] -> GetEntries());
00116 if(histEntry > 100){
00117 deltx = m_hxt[lay][iEntr][iLR][bin] -> GetMean();
00118 xerr = m_hxt[lay][iEntr][iLR][bin]->GetRMS();
00119 } else{
00120 continue;
00121 }
00122
00123 if(bin < NXTBIN)
00124 tbcen = ( (double)bin + 0.5 ) * gTbinw;
00125 else tbcen = xtini[6];
00126 xcor = xtFun(tbcen, xtini) - deltx;
00127
00128 if((tbcen <= Tmax) || (bin == NXTBIN)){
00129 TBINCEN.push_back( tbcen );
00130 XMEAS.push_back( xcor );
00131 ERR.push_back( xerr );
00132 } else{
00133 TBINCENED.push_back( tbcen );
00134 XMEASED.push_back( xcor );
00135 ERRED.push_back( xerr );
00136 }
00137 fxtlog << setw(3) << bin
00138 << setw(15) << deltx
00139 << setw(15) << xcor
00140 << setw(15) << tbcen
00141 << setw(15) << xerr
00142 << endl;
00143 }
00144
00145 if( XMEAS.size() < 12 ){
00146 TBINCEN.clear();
00147 XMEAS.clear();
00148 ERR.clear();
00149
00150 TBINCENED.clear();
00151 XMEASED.clear();
00152 ERRED.clear();
00153
00154 continue;
00155 }
00156
00157 for(int ord=0; ord<=5; ord++){
00158 arglist[0] = ord + 1;
00159 arglist[1] = xtini[ord];
00160 gmxt -> mnexcm("SET PARameter", arglist, 2, ierflg);
00161 }
00162
00163
00164 if(1 == gfixXtC0){
00165 arglist[0] = 1;
00166 arglist[1] = 0.0;
00167 gmxt -> mnexcm("SET PARameter", arglist, 2, ierflg);
00168 gmxt -> mnexcm("FIX", arglist, 1, ierflg);
00169 }
00170
00171 arglist[0] = 1000;
00172 arglist[1] = 0.1;
00173 gmxt -> mnexcm("MIGRAD", arglist, 2, ierflg);
00174 gmxt -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
00175
00176 fxtlog << "Xtpar: " << endl;
00177 if( (0 == ierflg) && (istat >= 2) ){
00178 for(int ord=0; ord<=5; ord++){
00179 gmxt -> GetParameter(ord, xtpar, xterr);
00180
00181 xtfit[ord] = xtpar;
00182
00183 if(1 == gNEntr[lay]){
00184 for(i=0; i<18; i++)
00185 calconst -> resetXtpar(lay, i, iLR, ord, xtpar);
00186 } else if(2 == gNEntr[lay]){
00187 if(0 == iEntr){
00188 for(i=0; i<9; i++)
00189 calconst->resetXtpar(lay, i, iLR, ord, xtpar);
00190 } else{
00191 for(i=9; i<18; i++)
00192 calconst->resetXtpar(lay, i, iLR, ord, xtpar);
00193 }
00194 }
00195 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
00196 }
00197 } else{
00198 for(int ord=0; ord<=5; ord++){
00199 fxtlog << setw(15) << xtini[ord] << setw(15) << "0" << endl;
00200 }
00201 }
00202 fxtlog << setw(15) << Tmax << setw(15) << "0" << endl;
00203
00204
00205 if(1 == gfixXtC0){
00206 arglist[0] = 1;
00207 gmxt -> mnexcm("REL", arglist, 1, ierflg);
00208 }
00209
00210 Dmax = xtFun(Tmax, xtfit);
00211
00212 if( XMEASED.size() >= 3 ){
00213
00214 arglist[0] = 1;
00215 arglist[1] = xtini[7];
00216 gmxtEd -> mnexcm("SET PARameter", arglist, 2, ierflg);
00217
00218 arglist[0] = 1000;
00219 arglist[1] = 0.1;
00220 gmxtEd -> mnexcm("MIGRAD", arglist, 2, ierflg);
00221 gmxtEd -> mnstat(fmin, edm, errdef, nvpar, nparx, istat);
00222
00223 if( (0 == ierflg) && (istat >=2) ){
00224 gmxtEd -> GetParameter(0, xtpar, xterr);
00225 if(xtpar < 0.0) xtpar = 0.0;
00226
00227
00228 if(1 == gNEntr[lay]){
00229 for(i=0; i<18; i++)
00230 calconst -> resetXtpar(lay, i, iLR, 7, xtpar);
00231 } else if(2 == gNEntr[lay]){
00232 if(0 == iEntr){
00233 for(i=0; i<9; i++)
00234 calconst->resetXtpar(lay, i, iLR, 7, xtpar);
00235 } else{
00236 for(i=9; i<18; i++)
00237 calconst->resetXtpar(lay, i, iLR, 7, xtpar);
00238 }
00239 }
00240 fxtlog << setw(15) << xtpar << setw(15) << xterr << endl;
00241 } else {
00242 fxtlog << setw(15) << xtini[7] << setw(15) << "0" << endl;
00243 }
00244 } else {
00245 fxtlog << setw(15) << xtini[7] << setw(15) << "0" << endl;
00246 }
00247 fxtlog << "Tm " << setw(15) << Tmax
00248 << " Dmax " << setw(15) << Dmax << endl;
00249
00250 TBINCEN.clear();
00251 XMEAS.clear();
00252 ERR.clear();
00253 TBINCENED.clear();
00254 XMEASED.clear();
00255 ERRED.clear();
00256 }
00257 }
00258 }
00259 fxtlog.close();
00260
00261 renameHist();
00262 delete gmxt;
00263 delete gmxtEd;
00264 }
00265
00266 void XtCalib::renameHist(){
00267 char hname[200];
00268 m_fdXt->SetName("fdXt");
00269 for(int lay=0; lay<43; lay++){
00270 for(int iEntr=0; iEntr<gNEntr[lay]; iEntr++){
00271 for(int lr=0; lr<NLR; lr++){
00272 for(int bin=0; bin<=NXTBIN; bin++){
00273 sprintf(hname, "Hxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, lr, bin);
00274 m_hxt[lay][iEntr][lr][bin]->SetName(hname);
00275 }
00276 }
00277 }
00278 }
00279 }