StRoot: StFtpcCalibMaker/StFtpcLaserCalib.cc Source File

00001 // $Id: StFtpcLaserCalib.cc,v 1.7 2008/05/15 22:39:47 jcs Exp $ 00002 // 00003 // $Log: StFtpcLaserCalib.cc,v $ 00004 // Revision 1.7 2008/05/15 22:39:47 jcs 00005 // re-activate helix fit 00006 // 00007 // Revision 1.6 2007/01/22 13:08:15 jcs 00008 // replace cout with LOG 00009 // 00010 // Revision 1.5 2007/01/19 08:22:59 jcs 00011 // replace true, false with kTRUE, kFALSE 00012 // 00013 // Revision 1.4 2006/06/19 12:37:58 jcs 00014 // Use negative value of zfieldcage for FTPC east 00015 // 00016 // Revision 1.3 2006/04/04 10:57:05 jcs 00017 // Fix memory leak 00018 // 00019 // Revision 1.2 2006/03/15 15:13:56 jcs 00020 // add lines for listing CVS update info 00021 // 00022 00023 #include "StFtpcLaserCalib.hh" 00024 #include "StFtpcTrackMaker/StFtpcPoint.hh" 00025 #include "StFtpcTrackMaker/StFtpcTrack.hh" 00026 #include "StFtpcTrackMaker/StFtpcVertex.hh" 00027 #include "TObjArray.h" 00028 00029 #define HELIX_FIT 00030 00031 const Float_t rad2grad=180/TMath::Pi(); 00032 00033 // define functions etc. used for fit added on 11/21/2003 00034 static Double_t funcxz(float x,float z,Double_t *par); 00035 static Double_t funcyz(float y,float z,Double_t *par); 00036 static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); 00037 00038 // 00039 00040 // im prinzip unschoen !!! Warning : defined but not used !!! 00041 static int nhits; // needed by the fcn used for minuit. 00042 static Float_t z[11],x[11],y[11]; 00043 static Float_t ex[11],ey[11]; 00044 00045 //______________________________________________________________________________ 00046 00047 StFtpcLaserCalib::StFtpcLaserCalib() 00048 { 00049 mMinuit=new TMinuit(4); 00050 //mMinuit->SetPrintLevel(-1); 00051 } 00052 00053 //______________________________________________________________________________ 00054 00055 //StFtpcLaserCalib::StFtpcLaserCalib(int ftpc, int lsec, int straight, int gfit, int minz, int maxz, int minrad, int maxrad, float gt0, float ggas, StFtpcLaserTrafo *trafo, StMagUtilities *gmagf) 00056 StFtpcLaserCalib::StFtpcLaserCalib(int ftpc, int lsec, int straight, int gfit, int minz, int maxz, int minrad, int maxrad, float gt0, float ggas, StFtpcLaserTrafo *trafo, StarMagField *gmagf) 00057 { 00058 FTPC=ftpc; 00059 LSEC=lsec; 00060 STRAIGHT=straight; 00061 MINZ=minz; 00062 MAXZ=maxz; 00063 MINRAD=minrad; 00064 MAXRAD=maxrad; 00065 GAUSFIT=gfit; 00066 mMinuit=new TMinuit(4); 00067 deltat0=gt0; 00068 deltagas=ggas; 00069 00070 usedfit=0; 00071 00072 #ifdef HELIX_FIT 00073 trcharge=0; 00074 p=pt=invp=invpt=0; 00075 #endif 00076 00077 magf=gmagf; 00078 00079 if (deltat0!=0 || deltagas!=0) 00080 mtrafo=trafo; 00081 } 00082 00083 //______________________________________________________________________________ 00084 00085 //static 00086 Double_t funcxz(float x,float z,Double_t *par) 00087 { 00088 00089 return x-(par[0]*z+par[1]); 00090 } 00091 00092 //static 00093 Double_t funcyz(float y,float z,Double_t *par) 00094 { 00095 return y-(par[2]*z+par[3]); 00096 } 00097 00098 //______________________________________________________________________________ 00099 00100 //static 00101 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) 00102 //void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) 00103 { 00104 00105 Int_t i; //calculate chisquare 00106 00107 Double_t chisq = 0; 00108 Double_t deltaxz, deltayz; 00109 //LOG_DEBUG<<"in fcn vor for .."<<nhits<<" hits"<<endm; 00110 00111 for (i=0;i<nhits; i++) 00112 { 00113 // Fehler aus Daten noch intgerieren ! 00114 // so richtig !?????? 00115 deltaxz=funcxz(x[i],z[i],par)/ex[i]; 00116 deltayz=funcyz(y[i],z[i],par)/ey[i]; 00117 chisq += (deltaxz*deltaxz+deltayz*deltayz); 00118 } 00119 00120 //LOG_DEBUG<<" f = "<<f<<endm; 00121 f = chisq; 00122 } 00123 00124 //______________________________________________________________________________ 00125 00126 bool StFtpcLaserCalib::cut(int i) 00127 { 00128 //LOG_DEBUG<<"in cut mit "<< z[i]<< " "<<hsec[i]<<endm; 00129 if ((z[i]>MINZ && z[i]<MAXZ) && (radius[i]>MINRAD && radius[i]<MAXRAD) && laser_sector(FTPC,LSEC,hsec[i]) ) 00130 return kTRUE; 00131 else 00132 return kFALSE; 00133 } 00134 00135 //______________________________________________________________________________ 00136 00137 void StFtpcLaserCalib::minuit_init() 00138 { 00139 //LOG_DEBUG<<"Minuit init ..."<<endm; 00140 //mMinuit = new TMinuit(4); //???? why is there a link error !? 00141 mMinuit->SetPrintLevel(-1); 00142 mMinuit->SetFCN(&fcn); 00143 //LOG_DEBUG<<"set fcn"<<endm; 00144 ierflg = 0; 00145 arglist[0] = 1; 00146 00147 mMinuit->mnexcm("SET ERR", arglist ,1,ierflg); 00148 //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm; 00149 } 00150 00151 //______________________________________________________________________________ 00152 00153 void StFtpcLaserCalib::minuit_set_par() 00154 { 00155 mMinuit->SetPrintLevel(-1); 00156 //LOG_DEBUG<<"Minuit set parameter ..."<<endm; 00157 //static Double_t vstart[4] = {1,1,1,1}; 00158 Double_t vstart[4] = {1,1,1,1}; 00159 // ggf. daten werte als Startparameter benutzen ! 00160 //static Double_t step[4] = {0.1,0.1,0.1,0.1}; 00161 Double_t step[4] = {0.1,0.1,0.1,0.1}; 00162 mMinuit->mnparm(0, "mxz", vstart[0], step[0], 0,0,ierflg); 00163 mMinuit->mnparm(1, "bxz", vstart[1], step[1], 0,0,ierflg); 00164 mMinuit->mnparm(2, "myz", vstart[2], step[2], 0,0,ierflg); 00165 mMinuit->mnparm(3, "byz", vstart[3], step[3], 0,0,ierflg); 00166 //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm; 00167 } 00168 00169 //______________________________________________________________________________ 00170 00171 void StFtpcLaserCalib::minuit_print() 00172 { 00173 // Print results 00174 Double_t amin,edm,errdef; 00175 Int_t nvpar,nparx,icstat; 00176 LOG_INFO<<"StFtpcLaserCalib::minuit_print():"<<endm; 00177 mMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat); 00178 mMinuit->mnprin(5,amin); 00179 LOG_INFO<<" "<<endm; 00180 } 00181 00182 //______________________________________________________________________________ 00183 00184 00185 void StFtpcLaserCalib::calc_res() 00186 { 00187 //LOG_DEBUG<<"Calculate residuals in x,y & r,phi ..."<<endm; 00188 //LOG_DEBUG<<" "<<endm; 00189 00190 00191 TF1 *fxz=new TF1("fxz","[0]*x+[1]",0,6); 00192 TF1 *fyz=new TF1("fyz","[0]*x+[1]",0,6); 00193 00194 // set fitted parameter to lin. function 00195 Double_t v,ev; 00196 for (int i=0;i<4;i++) 00197 { 00198 mMinuit->GetParameter(i,v,ev); 00199 if (i<2) fxz->SetParameter(i,v); 00200 else fyz->SetParameter(i-2,v); 00201 } 00202 00203 // calc residuals in x and y 00204 00205 hnhits->Fill(nhits); 00206 #ifdef HELIX_FIT 00207 htrcharge->Fill(trcharge); 00208 hp->Fill(p); 00209 hpt->Fill(pt); 00210 hinvp->Fill(1/p*trcharge); 00211 hinvpt->Fill(1/pt*trcharge); 00212 #endif 00213 00214 for (int i=0;i<nhits;i++) 00215 { 00216 calcx[i]=fxz->Eval(z[i]); 00217 calcy[i]=fyz->Eval(z[i]); 00218 00219 usedfit++; 00220 00221 resx[i]=(x[i]-fxz->Eval(z[i])); 00222 resy[i]=(y[i]-fyz->Eval(z[i])); 00223 calcrad[i]=sqrt(fxz->Eval(z[i])*fxz->Eval(z[i])+fyz->Eval(z[i])*fyz->Eval(z[i])); 00224 if (FTPC==1) 00225 rpol[i]=sqrt(fxz->Eval(zfieldcage)*fxz->Eval(zfieldcage)+fyz->Eval(zfieldcage)*fyz->Eval(zfieldcage)); 00226 if (FTPC==2) 00227 rpol[i]=sqrt(fxz->Eval(-zfieldcage)*fxz->Eval(-zfieldcage)+fyz->Eval(-zfieldcage)*fyz->Eval(-zfieldcage)); 00228 calcphi[i]=zyltrafo(fxz->Eval(z[i]),fyz->Eval(z[i]),z[i]); 00229 resrad[i]=radius[i]-calcrad[i]; 00230 resphi[i]=(phi[i]-calcphi[i])*Rad2Grad; 00231 00232 } 00233 00234 delete fxz; 00235 delete fyz; 00236 } 00237 00238 //______________________________________________________________________________ 00239 00240 00241 int StFtpcLaserCalib::laser_fit(int getnhits) 00242 { 00243 00244 nhits=getnhits; 00245 00246 int new_nhits=0; 00247 00248 for (int i=0;i<nhits;i++) 00249 { 00250 if (cut(i)) 00251 { 00252 fillarray(x[i],y[i],z[i],ex[i],ey[i],new_nhits,hsec[i],ppos[i],ppossigma[i],softsec[i],softrow[i],timepos[i],padl[i],timel[i],maxadc[i],charge[i]); 00253 new_nhits++; 00254 } 00255 } 00256 00257 nhits=new_nhits; 00258 if (nhits==0) return -1; 00259 00260 // *************************************************** 00261 // * calculate new space point depending on gas & t0 * 00262 // *************************************************** 00263 00264 if ((deltat0!=0 || deltagas!=0) && nhits>0) 00265 { 00266 //LOG_DEBUG<<"calculating new space point..."<<endm; 00267 for (int i=0;i<nhits;i++) // debug : <=nhits 00268 { 00269 00270 //LOG_DEBUG<<sqrt(x[i]*x[i]+y[i]*y[i])<<endm; 00271 //LOG_DEBUG<<softrow[i]<<" "<<softsec[i]<<" "<<timepos[i]<<" "<<ppos[i]<<" "<<x[i]<<" "<<y[i]<<endm; 00272 00273 if (mtrafo->padtrans(softrow[i],softsec[i],timepos[i],ppos[i],&x[i],&y[i])) 00274 { 00275 00276 radius[i]=sqrt(x[i]*x[i]+y[i]*y[i]); 00277 //LOG_DEBUG<<"radius["<<i<<"] = "<<radius[i]<<endm; 00278 phi[i]=zyltrafo(x[i],y[i],z[i]); 00279 } 00280 else { 00281 LOG_ERROR<<"ERROR padtrans !"<<endm; 00282 } 00283 } 00284 //LOG_DEBUG<<"after new space point calculation !"<<endm; 00285 } 00286 00287 #ifdef HELIX_FIT 00288 // ************************************************************************ 00289 // * include Helix fit & B-field distortion corrections (StFtpcTrack) !!!) * 00290 // * keine Trafo noetig da nur im Global, ggf. fuer East !??? * 00291 // ************************************************************************ 00292 00293 // create Laser Dummy Vertex (fuer data aus Tree auslesen !!! oder global nur); 00294 00295 StFtpcVertex *laser_vertex=new StFtpcVertex(0,0,0,0,0,0); 00296 00297 // create new track 00298 00299 StFtpcTrack *mLaserTrack=new StFtpcTrack(1); 00300 00301 TObjArray *mTrackHits=new TObjArray(0); 00302 00303 // create new point and add to track 00304 for (int i=0;i<nhits;i++) 00305 { 00306 // StFtpcPoint *mNewPoint=new StFtpcPoint(softrow[i],softsec[i],0,0,0,0,x[i],y[i],z[i],ex[i],ey[i],0,0,0,0); 00307 StFtpcPoint *mNewPoint=new StFtpcPoint(softrow[i],softsec[i],0,0,0,0,0,0,0,0,x[i],y[i],z[i],ex[i],ey[i],0,0,0,0); 00308 mTrackHits->AddLast(mNewPoint); 00309 //LOG_DEBUG<<mNewPoint->GetX()<<" "<<mNewPoint->GetY()<<" "<<mNewPoint->GetZ()<<endm; 00310 mLaserTrack->AddPoint((StFtpcPoint*) mTrackHits->Last()); 00311 //delete mNewPoint; 00312 } 00313 00314 // DEBUG : 00315 //LOG_DEBUG<<mLaserTrack->GetTrackNumber()<<" "<<mLaserTrack->GetNumberOfPoints()<<endm; 00316 //LOG_DEBUG<<laser_vertex->GetX()<<" "<<laser_vertex->GetY()<<" "<<laser_vertex->GetZ()<<endm; 00317 00318 // ******************************** 00319 // * Fit laser track with Helix ! * 00320 // ******************************** 00321 00322 // mLaserTrack->Fit(laser_vertex,10.,true); 00323 mLaserTrack->Fit(laser_vertex,10.,false); 00324 00325 // DEBUG : 00326 //LOG_DEBUG<<"p = "<<mLaserTrack->GetP()<<" | pt = "<<mLaserTrack->GetPt()<<" | charge = "<<mLaserTrack->GetCharge()<<endm; 00327 00328 // get momentum and charge sign of Helix Fit !!! 00329 trcharge=mLaserTrack->GetCharge(); 00330 p=mLaserTrack->GetP(); 00331 pt=mLaserTrack->GetPt(); 00332 //pt=1/(mLaserTrack->curvature()); 00333 //LOG_DEBUG<<mLaserTrack->curvature()<<endm; 00334 00335 // delete Hit Array and get residuals before ... 00336 for (int i=0;i<nhits;i++) 00337 { 00338 00339 resx2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetXGlobResidual(); 00340 resy2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetYGlobResidual(); 00341 resrad2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetRGlobResidual(); 00342 resphi2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetPhiGlobResidual(); 00343 00344 // DEBUG : 00345 //LOG_DEBUG<<" i = "<<i<<endm; 00346 //LOG_DEBUG<<((StFtpcPoint*) mTrackHits->At(i))->GetXS()<<" "<<((StFtpcPoint*) mTrackHits->At(i))->GetX()<<endm; 00347 //LOG_DEBUG<<((StFtpcPoint*) mTrackHits->At(i))->GetZS()<<" "<<((StFtpcPoint*) mTrackHits->At(i))->GetZ()<<endm; 00348 // noch nicht alle Punkte richtig !!!! 00349 // letzter Punkt wird nicht gesetzt (kleinstes z) vgl StFtpcTrack.cc !!! 00350 // nehme dann standard !!! 00351 00352 //LOG_DEBUG<<x[i]-((StFtpcPoint*) mTrackHits->At(i))->GetXS()<<endm; 00353 00354 // create shiftet space position !!! 00355 00356 // if (i<nhits-1) 00357 // { 00358 // x_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetXS(); 00359 // y_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetYS(); 00360 // } 00361 // else 00362 // { 00363 x_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetX(); 00364 y_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetY(); 00365 // } 00366 00367 //LOG_DEBUG<<x_s[i]-x[i]<<endm; 00368 00369 // calc difference as check !!! 00370 x_d[i]=x[i]-x_s[i];y_d[i]=y[i]-y_s[i]; 00371 00372 delete (StFtpcPoint*) mTrackHits->At(i); 00373 } 00374 00375 delete mTrackHits; 00376 00377 delete laser_vertex; 00378 delete mLaserTrack; 00379 00380 // *********************************** 00381 // * Straight line Fit with Miniut ! * 00382 // *********************************** 00383 00384 00385 // if gas or t0 changes take B-field corrected space pos. 00386 // (check if chain default = corrected !?) 00387 00388 //if ((deltat0!=0 || deltagas!=0) && nhits>0) 00389 00390 00391 for (int i=0;i<nhits;i++) 00392 { 00393 x[i]=x_s[i]; 00394 y[i]=y_s[i]; 00395 } 00396 00397 #endif 00398 00399 minuit_init(); 00400 minuit_set_par(); 00401 //mMinuit->SetPrintLevel(-1); 00402 00403 arglist[0] = 500; 00404 arglist[1] = 1; 00405 00406 if (nhits>0) 00407 { 00408 //LOG_DEBUG<<"start fit"<<endm; 00409 mMinuit->mnexcm("MIGRAD", arglist , 2 ,ierflg); 00410 //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm; 00411 //check if fit okay ! 00412 if (ierflg==0) 00413 { 00414 //minuit_print(); 00415 calc_res(); 00416 00417 // Fill histograms !!!! 00418 for (int i=0;i<nhits;i++) 00419 Fill(i); 00420 00421 } 00422 } 00423 else 00424 { 00425 LOG_ERROR<<"StFtpcLaserCalib::laser_fit - Number cluster on tracks = 0 !"<<endm; 00426 ierflg=-1; 00427 } 00428 00429 //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm; 00430 00431 nhits=0; 00432 return ierflg; 00433 } 00434 00435 //______________________________________________________________________________ 00436 00437 // erstmal test fuer Spectrum ! 00438 void StFtpcLaserCalib::defl_angle_transv() 00439 { 00440 TSpectrum *specrad=new TSpectrum(); 00441 //LOG_DEBUG<<"*******************************"<<endm; 00442 //LOG_DEBUG<<"Anzahl der Hits (rad) : "<<specrad->Search(hrad,3)<<endm; 00443 LOG_INFO<<"*******************************"<<endm; 00444 for (int i=0;i<specrad->Search(hrad,2.5);i++) 00445 { 00446 LOG_INFO<<"Radius "<<i<<" : "<<specrad->GetPositionX()[i]<<endm; 00447 } 00448 LOG_INFO<<"*******************************"<<endm; 00449 specrad->Delete(); 00450 } 00451 00452 //______________________________________________________________________________ 00453 00454 void StFtpcLaserCalib::defl_histograms_st() 00455 { 00456 char name[30]; 00457 char name2[30]; 00458 for (int i=1;i<11;i++) 00459 { 00460 sprintf(name,"rad_l%d",i); 00461 sprintf(name2,"pad_l%d",i); 00462 if (i<4) 00463 { 00464 hpadcut[i]=new TH1F(name2,name2,200,105,125); 00465 hradcut[i]=new TH1F(name,name,(int)((radcutmaxst[i]-radcutminst[i])/0.1),radcutminst[i],radcutmaxst[i]); 00466 } 00467 else 00468 { 00469 hpadcut[i]=new TH1F(); 00470 hradcut[i]=new TH1F(); 00471 } 00472 } 00473 } 00474 00475 //______________________________________________________________________________ 00476 00477 void StFtpcLaserCalib::defl_histograms() 00478 { 00479 char name[30]; 00480 char name2[30]; 00481 for (int i=1;i<11;i++) 00482 { 00483 sprintf(name,"rad_l%d",i); 00484 sprintf(name2,"pad_l%d",i); 00485 00486 hpadcut[i]=new TH1F(name2,name2,300,120,150); 00487 hradcut[i]=new TH1F(name,name,(int)((radcutmax[i]-radcutmin[i])/0.1),radcutmin[i],radcutmax[i]); 00488 } 00489 } 00490 00491 //______________________________________________________________________________ 00492 00493 00494 void StFtpcLaserCalib::extrapol_histograms() 00495 { 00496 char name[30]; 00497 00498 for (int i=1;i<4;i++) 00499 { 00500 sprintf(name,"radpol_l%d",i); 00501 hradpolcut[i]=new TH1F(name,name,(int)((radcutmaxst[i]-radcutminst[i])/0.025),radcutminst[i],radcutmaxst[i]); 00502 } 00503 } 00504 00505 //______________________________________________________________________________ 00506 00507 void StFtpcLaserCalib::fill_extrapol_histograms(float getradpol) 00508 { 00509 for (int i=1;i<4;i++) 00510 { 00511 if (getradpol>radcutminst[i] && getradpol<radcutmaxst[i]) 00512 { 00513 ((TH1F*) hradpolcut[i])->Fill(getradpol); 00514 } 00515 } 00516 } 00517 00518 //______________________________________________________________________________ 00519 00520 void StFtpcLaserCalib::fill_defl_histograms(float getrad, float getpadpos) 00521 { 00522 for (int i=1;i<11;i++) 00523 { 00524 if (getrad>radcutmin[i] && getrad<radcutmax[i]) 00525 { 00526 ((TH1F*) hradcut[i])->Fill(getrad); 00527 ((TH1F*) hpadcut[i])->Fill(getpadpos); 00528 00529 } 00530 } 00531 } 00532 00533 //______________________________________________________________________________ 00534 00535 00536 void StFtpcLaserCalib::fill_defl_histograms_st(float getrad, float getpadpos) 00537 { 00538 for (int i=1;i<4;i++) 00539 { 00540 if (getrad>radcutminst[i] && getrad<radcutmaxst[i]) 00541 { 00542 ((TH1F*) hradcut[i])->Fill(getrad); 00543 ((TH1F*) hpadcut[i])->Fill(getpadpos); 00544 00545 } 00546 } 00547 } 00548 00549 //______________________________________________________________________________ 00550 00551 void StFtpcLaserCalib::analyse_defl() 00552 { 00553 ofstream padfile,padfile2; 00554 padfile.open(filename+"_padpos.log",ios::out); 00555 padfile2.open(filename+"_padpos2.log",ios::out); 00556 00557 // use mean 00558 padfile<<"Padposition and radius with mean of histograms :"<<endl; 00559 padfile<<"------------------------------------------------"<<endl; 00560 padfile<<endl; 00561 padfile<<"radius [cm] | padposition"<<endl; 00562 padfile<<"-------------------------"<<endl; 00563 padfile<<endl; 00564 00565 for (int i=1;i<11;i++) 00566 { 00567 if (((TH1F*) hradcut[i])->GetEntries()>0) 00568 { 00569 padfile<<((TH1F*) hradcut[i])->GetMean()<<"+-"<<((TH1F*) hradcut[i])->GetRMS()<<" "; 00570 padfile<<((TH1F*) hpadcut[i])->GetMean()<<"+-"<<((TH1F*) hpadcut[i])->GetRMS()<<endl; 00571 padfile2<<((TH1F*) hradcut[i])->GetMean()<<" "<<((TH1F*) hradcut[i])->GetRMS()<<" "; 00572 padfile2<<((TH1F*) hpadcut[i])->GetMean()<<" "<<((TH1F*) hpadcut[i])->GetRMS()<<endl; 00573 } 00574 } 00575 00576 padfile<<endl; 00577 00578 padfile<<"Padposition and radius with mean of gausfit :"<<endl; 00579 padfile<<"------------------------------------------------"<<endl; 00580 padfile<<endl; 00581 padfile<<"radius [cm] | padposition"<<endl; 00582 padfile<<"-------------------------"<<endl; 00583 padfile<<endl; 00584 00585 // use gausfit !!! 00586 TF1 *gausfit3=new TF1("gausfit3","gaus"); 00587 00588 for (int i=1;i<11;i++) 00589 { 00590 if (((TH1F*) hradcut[i])->GetEntries()>0) 00591 { 00592 ((TH1F*) hradcut[i])->Fit("gausfit3","NQ"); 00593 padfile<<gausfit3->GetParameter(1)<<"+-"<<gausfit3->GetParameter(2)<<" "; 00594 padfile2<<gausfit3->GetParameter(1)<<" "<<gausfit3->GetParameter(2)<<" "; 00595 ((TH1F*) hpadcut[i])->Fit("gausfit3","NQ"); 00596 padfile<<gausfit3->GetParameter(1)<<"+-"<<gausfit3->GetParameter(2)<<endl; 00597 padfile2<<gausfit3->GetParameter(1)<<" "<<gausfit3->GetParameter(2)<<endl; 00598 } 00599 } 00600 00601 // make ps file !!! 00602 00603 00604 padfile.close(); 00605 padfile2.close(); 00606 } 00607 00608 //______________________________________________________________________________ 00609 00610 void StFtpcLaserCalib::PositionLog() 00611 { 00612 00613 ofstream logfile; 00614 logfile.open(filename+".log",ios::out); 00615 00616 TSpectrum *spec=new TSpectrum(); 00617 //LOG_DEBUG<<" "<<endm; 00618 //LOG_DEBUG<<"Gaus fit to radius distribution of the "<<spec->Search(hrad,6)<<" reconstructed straight Laser tracks..."<<endm; 00619 //LOG_DEBUG<<"-----------------------------------------------------------------------------"<<endm; 00620 for (int i=0;i<spec->Search(hrad,6);i++) 00621 { 00622 //LOG_DEBUG<<"Radius after Peakfinder : "<<(spec->GetPositionX())[i]<<endm; 00623 logfile<<"Radius after Peakfinder : "<<(spec->GetPositionX())[i]<<endl; 00624 if (GAUSFIT==1) 00625 { 00626 float min=(spec->GetPositionX()[i]-1); 00627 float max=(spec->GetPositionX()[i]+1); 00628 TF1 *gausfit=new TF1("gausfit","gaus",min,max); 00629 gausfit->SetParameter(1,(spec->GetPositionX())[i]); 00630 // ggf. ParLimits einfuehren um richtige Fehlerbehandlung !!! 00631 hrad->Fit(gausfit,"RQN"); 00632 //LOG_DEBUG<<"Radius after Gausfit : "<<gausfit->GetParameter(1)<<" +- "<<gausfit->GetParameter(2)<<endm; 00633 logfile<<"Radius after Gausfit : "<<gausfit->GetParameter(1)<<" +- "<<gausfit->GetParameter(2)<<endl; 00634 } 00635 } 00636 //LOG_DEBUG<<" "<<endm; 00637 //LOG_DEBUG<<"Gaus fit to phi distribution of the "<<spec->Search(hphi,2)<<" reconstructed straight Laser tracks..."<<endm; 00638 //LOG_DEBUG<<"-----------------------------------------------------------------------------"<<endm; 00639 for (int i=0;i<spec->Search(hphi,6);i++) 00640 { 00641 //LOG_DEBUG<<"Phi after Peakfinder : "<<(spec->GetPositionX())[i]<<endm; 00642 logfile<<"Phi after Peakfinder : "<<(spec->GetPositionX())[i]<<endl; 00643 if (GAUSFIT==1) 00644 { 00645 float min=(spec->GetPositionX()[i]-0.5); 00646 float max=(spec->GetPositionX()[i]+0.5); 00647 TF1 *gausfit2=new TF1("gausfit2","gaus",min,max); 00648 gausfit2->SetParameter(1,(spec->GetPositionX())[i]); 00649 hphi->Fit(gausfit2,"RQN"); 00650 //LOG_DEBUG<<"Phi after Gausfit : "<<gausfit2->GetParameter(1)<<" +- "<<gausfit2->GetParameter(2)<<endm; 00651 logfile<<"Phi after Gausfit : "<<gausfit2->GetParameter(1)<<" +- "<<gausfit2->GetParameter(2)<<endl; 00652 } 00653 } 00654 logfile.close(); 00655 } 00656 00657 //______________________________________________________________________________ 00658 00659 void StFtpcLaserCalib::fillarray(float tx,float ty,float tz,float tex,float tey,int n,int nsec,float gppos,float gppossigma, int gsoftsec, int gsoftrow, float gtimepos, float getpadl, float gettimel, float getmaxadc, float getcharge) 00660 { 00661 radius[n]=sqrt(tx*tx+ty*ty); 00662 phi[n]=zyltrafo(tx,ty,tz); 00663 //phi[n]=atan2(ty,tx); 00664 x[n]=tx; 00665 y[n]=ty; 00666 z[n]=tz; 00667 hsec[n]=nsec; 00668 ex[n]=tex; 00669 ey[n]=tey; 00670 ppos[n]=gppos; 00671 ppossigma[n]=gppossigma; 00672 softsec[n]=gsoftsec; 00673 softrow[n]=gsoftrow; 00674 timepos[n]=gtimepos; 00675 timel[n]=gettimel; 00676 padl[n]=getpadl; 00677 charge[n]=getcharge; 00678 maxadc[n]=getmaxadc; 00679 00680 //LOG_DEBUG<<radius[n]<<" "<<maxadc[n]<<endm; 00681 } 00682 00683 //______________________________________________________________________________ 00684 00685 void StFtpcLaserCalib::Fill(int l) 00686 { 00687 hrad->Fill(radius[l]); 00688 hradpol->Fill(rpol[l]); 00689 hphi->Fill(phi[l]*rad2grad); 00690 hcalcrad->Fill(calcrad[l]); 00691 hcalcphi->Fill(calcphi[l]); 00692 hphiz->Fill(z[l],phi[l]); 00693 hradz->Fill(z[l],radius[l]); 00694 hresx->Fill(resx[l]); 00695 hresy->Fill(resy[l]); 00696 hresrad->Fill(resrad[l]); 00697 hresphi->Fill(resphi[l]); 00698 hresrad2->Fill(resrad[l],radius[l]); 00699 hresphi2->Fill(resphi[l],radius[l]); 00700 hresx2->Fill(resx[l],radius[l]); 00701 hresy2->Fill(resy[l],radius[l]); 00702 00703 #ifdef HELIX_FIT 00704 hhresx->Fill(resx2[l]); 00705 hhresy->Fill(resy2[l]); 00706 hhresrad->Fill(resrad2[l]); 00707 hhresphi->Fill(resphi2[l]*rad2grad); 00708 hhresrad2->Fill(resrad2[l],radius[l]); // radius ggf. da Helix mit schieben anders !???? 00709 hhresphi2->Fill(resphi2[l]*rad2grad,radius[l]); 00710 hhresx2->Fill(resx2[l],radius[l]); 00711 hhresy2->Fill(resy2[l],radius[l]); 00712 #endif 00713 00714 hpad->Fill(ppos[l]); 00715 htime->Fill(timepos[l]); 00716 hpadrad->Fill(radius[l],ppos[l]); 00717 hpadsigma->Fill(ppossigma[l]); 00718 hpadl->Fill(padl[l]); 00719 htimel->Fill(timel[l]); 00720 hcharge->Fill(charge[l]); 00721 hmaxadc->Fill(maxadc[l]); 00722 00723 #ifdef HELIX_FIT 00724 hbdiffx->Fill(x_d[l]); 00725 hbdiffy->Fill(y_d[l]); 00726 #endif 00727 00728 if (STRAIGHT==0 || STRAIGHT==3) 00729 fill_defl_histograms(radius[l],ppos[l]); 00730 else if (STRAIGHT==1) 00731 fill_defl_histograms_st(radius[l],ppos[l]); 00732 00733 fill_extrapol_histograms(rpol[l]); 00734 } 00735 00736 00737 //______________________________________________________________________________ 00738 00739 void StFtpcLaserCalib::MakePs() 00740 { 00741 LOG_INFO<<"Create "<<filename+".ps"<<endm; 00742 00743 gStyle->SetPalette(1); 00744 00745 TCanvas *c1 = new TCanvas("c1","ps",200,10,700,500); 00746 TPostScript *fps=new TPostScript(filename+".ps",112); 00747 00748 00749 fps->NewPage(); 00750 c1->Divide(2,2); 00751 c1->Update(); 00752 c1->cd(1); 00753 hrad->DrawCopy(); 00754 c1->cd(2); 00755 //hphi->DrawCopy(); 00756 hnhits->DrawCopy(); 00757 c1->cd(3); 00758 hradz->DrawCopy("box"); 00759 c1->cd(4); 00760 hphi->DrawCopy(); 00761 //hphiz->DrawCopy("box"); 00762 c1->Update(); 00763 fps->NewPage(); 00764 c1->cd(1); 00765 00766 TF1 *gfit=new TF1("gfit","gaus",-0.5,0.5); 00767 gfit->FixParameter(1,0.0); 00768 //gfit->SetParameter(1,0.0); 00769 //gfit->SetParLimits(1,0.0,0.0); 00770 gfit->SetLineColor(2); 00771 gfit->SetLineWidth(1); 00772 00773 ofstream logfile; 00774 logfile.open(filename+"_res.log",ios::out); 00775 c1->cd(1); 00776 hresx->DrawCopy(); 00777 hresx->Fit(gfit,"QR"); 00778 logfile<<deltat0<<" "<<deltagas<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl; 00779 c1->cd(2); 00780 hresy->DrawCopy(); 00781 hresy->Fit(gfit,"QR"); 00782 logfile<<deltat0<<" "<<deltagas<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl; 00783 c1->cd(3); 00784 hresrad->DrawCopy(); 00785 hresrad->Fit(gfit,"QR"); 00786 logfile<<deltat0<<" "<<deltagas<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl; 00787 c1->cd(4); 00788 hresphi->DrawCopy(); 00789 hresphi->Fit(gfit,"QR"); 00790 logfile<<deltat0<<" "<<deltagas<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl; 00791 //logfile<<"Delta t0 | Delta Gas | Sigma res (x,y,r,phi)"<<endl; 00792 //logfile.close(); 00793 00794 c1->Update(); 00795 fps->NewPage(); 00796 c1->cd(1); 00797 hresrad2->DrawCopy("colz"); 00798 c1->cd(2); 00799 hresphi2->DrawCopy("colz"); 00800 //c1->cd(3); 00801 //hrad->DrawCopy();hcalcrad->DrawCopy("same"); 00802 //c1->cd(4); 00803 //hphi->DrawCopy();hcalcphi->DrawCopy("same"); 00804 c1->cd(3); 00805 hresx2->DrawCopy("colz"); 00806 c1->cd(4); 00807 hresy2->DrawCopy("colz"); 00808 c1->Update(); 00809 fps->NewPage(); 00810 c1->cd(1); 00811 hradpol->Draw(); 00812 //logfile<<deltat0<<" "<<deltagas<<" "<<hradpol->GetMean()<<endl; 00813 //logfile.close(); 00814 c1->cd(2); 00815 hpad->DrawCopy(); 00816 c1->cd(3); 00817 hpadrad->DrawCopy("colz"); 00818 c1->cd(4); 00819 hpadsigma->DrawCopy(); 00820 c1->Update(); 00821 fps->NewPage(); 00822 00823 for (int i=1;i<4;i++) 00824 { 00825 c1->cd(i); 00826 ((TH1F*) hradpolcut[i])->DrawCopy(); 00827 logfile<<deltat0<<" "<<deltagas<<" "<<((TH1F*) hradpolcut[i])->GetMean()<<" "<<((TH1F*) hradpolcut[i])->GetRMS()<<endl; 00828 logfile<<deltat0<<" "<<deltagas<<" "<<(((TH1F*) hradpolcut[i])->GetXaxis())->GetBinCenter(((TH1F*) hradpolcut[i])->GetMaximumBin())<<endl; 00829 } 00830 c1->Update(); 00831 00832 fps->NewPage(); 00833 c1->cd(1);hpadl->DrawCopy(); 00834 c1->cd(2);htimel->DrawCopy(); 00835 c1->cd(3);hmaxadc->DrawCopy();//hmaxadc->Fit("landau"); 00836 c1->cd(4);hcharge->DrawCopy();//hcharge->Fit("landau"); 00837 c1->Update(); 00838 00839 #ifdef HELIX_FIT 00840 fps->NewPage(); 00841 //c1->cd(1);htrcharge->DrawCopy(); 00842 c1->cd(1);hinvp->DrawCopy();hinvp->Fit(gfit,"QR"); 00843 c1->cd(2);hinvpt->DrawCopy();hinvpt->Fit(gfit,"QR"); 00844 c1->cd(3);hbdiffx->DrawCopy(); 00845 c1->cd(4);hbdiffy->DrawCopy(); 00846 c1->Update(); 00847 00848 gfit->SetLineColor(3); 00849 00850 fps->NewPage(); 00851 c1->cd(1); 00852 hhresx->DrawCopy(); 00853 hhresx->Fit(gfit,"QR"); 00854 c1->cd(2); 00855 hhresy->DrawCopy(); 00856 hhresy->Fit(gfit,"QR"); 00857 c1->cd(3); 00858 hhresrad->DrawCopy(); 00859 hhresrad->Fit(gfit,"QR"); 00860 c1->cd(4); 00861 hhresphi->DrawCopy(); 00862 hhresphi->Fit(gfit,"QR"); 00863 c1->Update(); 00864 00865 fps->NewPage(); 00866 c1->cd(1); 00867 hhresrad2->DrawCopy("colz"); 00868 c1->cd(2); 00869 hhresphi2->DrawCopy("colz"); 00870 c1->cd(3); 00871 hhresx2->DrawCopy("colz"); 00872 c1->cd(4); 00873 hhresy2->DrawCopy("colz"); 00874 c1->Update(); 00875 00876 fps->NewPage(); 00877 c1->cd(1); 00878 hhresx->DrawCopy(); 00879 hresx->SetLineColor(2); 00880 hresx->DrawCopy("same"); 00881 c1->cd(2); 00882 hhresy->DrawCopy(); 00883 hresy->SetLineColor(2); 00884 hresy->DrawCopy("same"); 00885 c1->cd(3); 00886 hhresrad->DrawCopy(); 00887 hresrad->SetLineColor(2); 00888 hresrad->DrawCopy("same"); 00889 c1->cd(4); 00890 hhresphi->DrawCopy(); 00891 hresphi->SetLineColor(2); 00892 hresphi->DrawCopy("same"); 00893 c1->Update(); 00894 #endif 00895 00896 logfile.close(); 00897 fps->Close(); 00898 //c1->Delete(); 00899 } 00900 00901 //______________________________________________________________________________ 00902 00903 void StFtpcLaserCalib::MakeOutput(TString eingabe,char* t0, char* gas) 00904 { 00905 00906 00907 eingabe +="-pos_"; 00908 if (FTPC==1) 00909 eingabe +="w_lsec_"; 00910 else if (FTPC==2) 00911 eingabe +="e_lsec_"; 00912 else if (FTPC==0) 00913 eingabe +="all_lsec"; 00914 00915 if (FTPC!=0) 00916 eingabe +=LSEC; 00917 00918 if (STRAIGHT==1) 00919 eingabe+="_g"; 00920 else if (STRAIGHT==0) 00921 eingabe+="_s"; 00922 else if (STRAIGHT==2) 00923 eingabe+="_s2"; 00924 else if (STRAIGHT==3) 00925 eingabe+="_all"; 00926 00927 eingabe +="_dt"; eingabe +=t0; 00928 eingabe +="_dg"; eingabe +=gas; 00929 00930 filename=eingabe; 00931 00932 LOG_INFO<<"Create file : "<<filename+".root"<<endm; 00933 outfile=new TFile(filename+".root","RECREATE"); 00934 00935 hrad=new TH1F("rad","radius (straight) laser tracks",124*8,0.5,31.5); 00936 hradpol=new TH1F("radpol","radius laser tracks extrapoliert fieldcage",124*8,0.5,31.5); 00937 //hphi=new TH1F("phi","phi (straight) laser tracks",48*8,-6,6); 00938 hphi=new TH1F("phi","phi laser tracks",360*4,-90,270); 00939 hcalcrad=new TH1F("clacrad","calc radius",124*8,0.5,31.5); 00940 hcalcphi=new TH1F("calcphi","calc phi",48*8,-6,6); 00941 hcalcrad->SetLineColor(2);hcalcphi->SetLineColor(2); 00942 00943 if (FTPC==1) 00944 { 00945 hradz=new TH2F("radz","radius (straight) laser tracks",100,140,275,124*4,0.5,31.5); 00946 hphiz=new TH2F("phiz","phi (straight) laser tracks",100,140,275,48*2,0,4); 00947 } 00948 else if (FTPC==2) 00949 { 00950 hradz=new TH2F("radz","radius (straight) laser tracks",100,-275,-140,124*4,0.5,31.5); 00951 hphiz=new TH2F("phiz","phi (straight) laser tracks",100,-275,-140,48*2,0,4); 00952 } 00953 else if (FTPC==0) 00954 { 00955 hradz=new TH2F("radz","radius (straight) laser tracks",200,-275,275,124*4,0.5,31.5); 00956 hphiz=new TH2F("phiz","phi (straight) laser tracks",200,-275,275,48*2,0,4); 00957 } 00958 00959 Int_t r_vs_z_bin=124/2; 00960 00961 hresx=new TH1F("resx","Residual x",100,-0.5,0.5); 00962 hresy=new TH1F("resy","Residual y",100,-0.5,0.5); 00963 hresrad=new TH1F("resrad","Residual radius",100,-0.5,0.5); 00964 hresphi=new TH1F("resphi","Residual phi",100,-0.5,0.5); 00965 hresx2=new TH2F("resx2","Residual x vs.radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 00966 hresx2->SetMinimum(0); 00967 hresy2=new TH2F("resy2","Residual y vs. radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 00968 hresy2->SetMinimum(0); 00969 hresrad2=new TH2F("resrad2","Residual radius vs.radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 00970 hresrad2->SetMinimum(0); 00971 hresphi2=new TH2F("resphi2","Residual phi vs. radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 00972 hresphi2->SetMinimum(0); 00973 hpad=new TH1F("pados","Padposition",1600,0,160); 00974 htime=new TH1F("timepos","Timeposition",1800,0,180); 00975 hpadrad=new TH2F("radpados","Padposition vs. radius",r_vs_z_bin,0.5,31.5,1600,0,160); 00976 hpadsigma=new TH1F("padsigma","Sigma padposition",50,0,5); 00977 00978 hmaxadc=new TH1F("maxadc","MaxAdc",150,0.5,150.5); 00979 hcharge=new TH1F("charge","Charge",300,0.5,1500.5); 00980 hpadl=new TH1F("padl","Padlength",15,0.5,15.5); 00981 htimel=new TH1F("timel","Timelength",20,0.5,20.5); 00982 00983 hnhits=new TH1F("nhits","Number hits on laser track",8,3.5,11.5); 00984 00985 #ifdef HELIX_FIT 00986 htrcharge=new TH1F("trcharge","Charge of Helix Fit",3,-1.5,1.5); 00987 hp=new TH1F("p","p Helix Fit",100, -0.5, 0.5); 00988 hpt=new TH1F("pt","pt Helix Fit",100, -0.5, 0.5); 00989 hinvp=new TH1F("invp","1/p * charge of Helix Fit",100, -0.5, 0.5); 00990 hinvpt=new TH1F("invpt","1/pt * charge of Helix Fit",50, -0.5, 0.5); 00991 00992 hhresx=new TH1F("hresx","Residual x Helix Fit",100,-0.5,0.5); 00993 hhresy=new TH1F("hrresy","Residual y Helix Fit",100,-0.5,0.5); 00994 hhresrad=new TH1F("hrresrad","Residual radius Helix Fit",100,-0.5,0.5); 00995 hhresphi=new TH1F("hresphi","Residual phi Helix Fit",100,-0.5,0.5); 00996 hhresrad2=new TH2F("hresrad2","Residual radius vs.radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 00997 hhresphi2=new TH2F("hresphi2","Residual phi vs. radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 00998 hhresx2=new TH2F("hresx2","Residual x vs.radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 00999 hhresy2=new TH2F("hresy2","Residual y vs. radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5); 01000 01001 hbdiffx=new TH1F("bdiffx","Differenz x and B-field corr. x",100,-0.01,0.01); 01002 hbdiffy=new TH1F("bdiffy","Differenz y and B-field corr. y",100,-0.01,0.01); 01003 // if you activate the following two histograms be sure to activate the corresponding delete statement 01004 //hbdiffx2=new TH2F("bdiffx2","Differenz x and B-field corr. x vs. z",); 01005 //hbdiffy2=new TH2F("bdiffy2","Differenz y and B-field corr. y vs. z",); 01006 #endif 01007 01008 if (STRAIGHT==0) 01009 defl_histograms(); 01010 else 01011 defl_histograms_st(); 01012 01013 //LOG_DEBUG<<"For extrapolhistos"<<endm; 01014 extrapol_histograms(); 01015 01016 } 01017 01018 01019 //______________________________________________________________________________ 01020 01021 StFtpcLaserCalib::~StFtpcLaserCalib() 01022 { 01023 //if (deltat0=0) 01024 //delete trafo; 01025 01026 outfile->Write(); 01027 01028 delete mMinuit; 01029 01030 //LOG_DEBUG<<"for normal histos"<<endm; 01031 01032 delete hrad;delete hradpol;delete hphi; delete hcalcrad; 01033 delete hcalcphi; delete hradz; delete hphiz; 01034 delete hresx;delete hresy;delete hresphi; delete hresrad; delete hresrad2; 01035 delete hresphi2; delete hpad; delete hpadrad; delete hpadsigma; delete htime; 01036 delete hnhits; 01037 #ifdef HELIX_FIT 01038 delete hhresx;delete hhresy;delete hhresphi; delete hhresrad; delete hhresrad2; delete hhresphi2; 01039 delete htrcharge; 01040 delete hp; delete hpt; delete hinvp; delete hinvpt; 01041 delete hhresx2; delete hhresy2; 01042 delete hresx2; delete hresy2; 01043 delete hbdiffx; delete hbdiffy; 01044 //delete hbdiffx2;delete hbdiffy2; 01045 #endif 01046 01047 if (STRAIGHT==0) 01048 { 01049 for (int i=1;i<4;i++) 01050 {delete hpadcut[i];delete hradcut[i];} 01051 } 01052 else 01053 { 01054 for (int i=1;i<11;i++) 01055 {delete hpadcut[i];delete hradcut[i];} 01056 } 01057 outfile->Close(); 01058 //LOG_DEBUG<<"StFtpcLaserCalib() deconstructed !"<<endm; 01059 }