StRoot: StGenericVertexMaker/macros/mainFitBeamLine3D.cxx Source File

00001 #include <stdio.h>
00002 #include <math.h>
00003 #include <assert.h>
00004 #include <TMinuit.h>
00005 #include <TFile.h>
00006 #include <TH1.h>
00007 #include <TObjArray.h>
00008 #include <TRandom3.h>
00009 
00010 #include "UtilBeamLine3D.h"
00011 
00012 //#define FIT_XYSLOPES
00013 // Do not use it for now, fit stops at strange minima
00014 
00015 UtilBeamLine3D util; // global instance
00016 int twoLineDca3D(double& lambda,double& kappa,TVector3& V, TVector3& U, TVector3& R,TVector3& P);
00017 
00018 //===========================
00019 int main() {
00020   printf("THIS is main\n");
00021   int mode = 3; //0 = just migrad, no tilt
00022                 //1 = migrad - position + tilt, but not both
00023                 //2 = migrad - position + tilt together
00024                 //3 = simplex-migrad-simplex-migrad
00025 
00026   const char *chararray[40];
00027   chararray[0] = "F10383";
00028   chararray[1] = "F10398";
00029   chararray[2] = "F10399";
00030   chararray[3] = "F10402";
00031   chararray[4] = "F10403";
00032   chararray[5] = "F10404";
00033   chararray[6] = "F10407";
00034   chararray[7] = "F10412";
00035   chararray[8] = "F10415";
00036   chararray[9] = "F10426";
00037   chararray[10] = "F10434";
00038   chararray[11] = "F10439";
00039   chararray[12] = "F10448";
00040   chararray[13] = "F10449";
00041   chararray[14] = "F10450";
00042   chararray[15] = "F10454";
00043   chararray[16] = "F10455";
00044   chararray[17] = "F10463";
00045   chararray[18] = "F10464";
00046   chararray[19] = "F10465";
00047   chararray[20] = "F10471";
00048   chararray[21] = "F10476";
00049   chararray[22] = "F10478";
00050   chararray[23] = "F10482";
00051   chararray[24] = "F10486";
00052   chararray[25] = "F10490";
00053   chararray[26] = "F10494";
00054   chararray[27] = "F10505";
00055   chararray[28] = "F10507";
00056   chararray[29] = "F10508";
00057   chararray[30] = "F10517";
00058   chararray[31] = "F10525";
00059   chararray[32] = "F10526";
00060   chararray[33] = "F10527";
00061   chararray[34] = "F10528";
00062   chararray[35] = "F10531";
00063   chararray[36] = "F10532";
00064   chararray[37] = "F10535";
00065   chararray[38] = "F10536";
00066   
00067 //   char *chararray[11];
00068 //   chararray[0] ="MC-pQCD-setE1"; // asorted M-C
00069 //   chararray[1] ="MC-pQCD-setE2"; // asorted M-C
00070 //   chararray[2] ="MC-pQCD-setE3"; // asorted M-C
00071 //   chararray[3] ="MC-pQCD-setE4"; // asorted M-C
00072 //   chararray[4] ="MC-pQCD-setE5"; // asorted M-C
00073 //   chararray[5] ="MC-pQCD-setE6"; // asorted M-C
00074 //   chararray[6] ="MC-pQCD-setE3+4"; // asorted M-C
00075 //   chararray[7] ="MC-pQCD-setE3+5"; // asorted M-C
00076 //   chararray[8] ="MC-pQCD-setE3+6"; // asorted M-C
00077 //   chararray[9] ="MC-pQCD-setE4+5"; // asorted M-C
00078 //   chararray[10] ="MC-pQCD-setE5+6"; // asorted M-C
00079   
00080   TRandom3* rnd = new TRandom3(1997);
00081   for (int i = 34;i<35;++i){
00082     
00083     //define variables needed to manipulate parameters after fit
00084     Double_t amin,edm,errdef;
00085     Int_t nvpar,nparx,icstat;
00086     double pval[mxPar],perr[mxPar],plo[mxPar],phi[mxPar];
00087     TString paraName[mxPar];
00088     int istat;
00089     const char *core = chararray[i];
00090     double dmax2=util.cut_Dmax*util.cut_Dmax; // to speed up 
00091 
00092     TObjArray  HList;
00093     util.initHisto(&HList);
00094 
00095     for (int j = 0;j<1;++j){ 
00096       pval[0] = (rnd->Rndm()-0.5)*2;
00097       pval[1] = (rnd->Rndm()-0.5)*2;
00098       pval[2] = (rnd->Rndm()-0.5)/25;
00099       pval[3] = (rnd->Rndm()-0.5)/25;
00100       //pval[2] = 0;
00101       //pval[3] = 0;
00102       
00103       printf("rand # starting point = %f , %f , %f , %f \n",pval[0],pval[1],pval[2],pval[3]);
00104       printf("main fit 3D beamline to %s start...\n",core);
00105       // select input file:
00106       TString inpFile="data/globTr_"; inpFile+=core; inpFile+=".rosi";
00107       util.readTracks(inpFile);
00108       
00109       util.print();
00110       util.qaTracks();
00111       
00112       TMinuit *gMinuit = new TMinuit(4);
00113       //..... pass name of the minimized function ......
00114       gMinuit->SetFCN(beamLineLike3D);
00115       double arglist[10];
00116       int error_flag = 0;
00117       
00118       arglist[0] = 1;
00119       gMinuit->mnexcm("SET ERR",arglist,1,error_flag);
00120       //Interprets a command and takes appropriate action
00121       //void mnexcm(const char* comand, Double_t* plist, Int_t llist, Int_t& ierflg)
00122       //double f_null = 0.0;
00123       
00124       //define range, step size and name of all parameters
00125       gMinuit->mnparm(0, "X0 (cm)",pval[0],0.01,-2,2,error_flag);
00126       gMinuit->mnparm(1, "Y0 (cm)",pval[1],0.01,-2,2,error_flag);
00127       gMinuit->mnparm(2, "Ux",0.0,0.001,-0.1,0.1,error_flag);
00128       gMinuit->mnparm(3, "Uy",0.0,0.001,-0.1,0.1,error_flag);
00129       
00130       
00131       arglist[0] = 500.;
00132       arglist[1] = 1.;
00133       
00134       //pick MINUIT algorithim to use
00135       
00136       const char* com = "MIGRAD";
00137       const char* com1 = "SIMPLEX";
00138       const char* com2 = "MIGRAD";
00139       //const char* com = "MINOS";
00140       
00141       /* ******* 
00142          minimize
00143          ******* */
00144       //fix tilt
00145       gMinuit->FixParameter(2);
00146       gMinuit->FixParameter(3);
00147       
00148       if (mode<3)
00149         gMinuit->mnexcm(com,arglist,2,error_flag);
00150       else
00151         {
00152           gMinuit->mnexcm(com1,arglist,2,error_flag);
00153           gMinuit->mnexcm(com2,arglist,2,error_flag);
00154           gMinuit->mnexcm(com1,arglist,2,error_flag);
00155           gMinuit->mnexcm(com2,arglist,2,error_flag);
00156         }
00157       
00158       gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
00159       gMinuit->mnprin(3,amin);
00160       //Prints the values of the parameters at the time of the call*-
00161       // 3    values, errors, step sizes, first derivs.
00162       
00163       util.fcnCount=-1; // disable ceratin monitorings for final plots
00164       
00165       
00166       for(int ip=0;ip<mxPar;ip++) {
00167         gMinuit->mnpout(ip,paraName[ip],pval[ip],perr[ip],plo[ip],phi[ip],istat);
00168         util.hA[19]->Fill(paraName[ip],pval[ip]);
00169         util.hA[19]->SetBinError(ip+1,perr[ip]);  
00170       }
00171       //printf("PARAValues X,Y,Ux,Uy for %s = %f , %f , %f , %f \n",core,pval[0],pval[1],pval[2],pval[3]);
00172       printf("PARAValues X,Y for %s = %f , %f \n",core,pval[0],pval[1]);
00173       if (1){
00174         if (mode > 0){
00175           printf("\n Now fitting tilt");
00176           
00177           TMinuit *tMinuit = new TMinuit(4);
00178           tMinuit->SetFCN(beamLineLike3D);
00179           arglist[0] = 1;
00180           tMinuit->mnexcm("SET ERR",arglist,1,error_flag);
00181           tMinuit->mnparm(0, "X0 (cm)",pval[0],0.01,-2,2,error_flag);
00182           tMinuit->mnparm(1, "Y0 (cm)",pval[1],0.01,-2,2,error_flag);
00183           tMinuit->mnparm(2, "Ux",pval[2],0.001,-0.02,0.02,error_flag);
00184           tMinuit->mnparm(3, "Uy",pval[3],0.001,-0.02,0.02,error_flag);
00185           
00186           arglist[0] = 500.;
00187           arglist[1] = 1.;
00188           
00189           tMinuit->FixParameter(0);
00190           tMinuit->FixParameter(1);
00191           
00192           if (mode < 3)
00193             tMinuit->mnexcm(com,arglist,2,error_flag);
00194           else
00195             {
00196               tMinuit->mnexcm(com1,arglist,2,error_flag);
00197               tMinuit->mnexcm(com2,arglist,2,error_flag);
00198               tMinuit->mnexcm(com1,arglist,2,error_flag);
00199               tMinuit->mnexcm(com2,arglist,2,error_flag);
00200             }
00201           tMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
00202           tMinuit->mnprin(3,amin);
00203           for(int ip=0;ip<mxPar;ip++) {
00204             tMinuit->mnpout(ip,paraName[ip],pval[ip],perr[ip],plo[ip],phi[ip],istat);}
00205           printf("PARAValues UX,UY for %s = %f , %f \n",core,pval[2],pval[3]);
00206           tMinuit->Delete();
00207         }
00208         if (mode > 1){
00209           printf("\n Now fitting both position and tilt");
00210           TMinuit *aMinuit = new TMinuit(4);
00211           aMinuit->SetFCN(beamLineLike3D);
00212           arglist[0] = 1;
00213           aMinuit->mnexcm("SET ERR",arglist,1,error_flag);
00214           aMinuit->mnexcm("SET ERR",arglist,1,error_flag);
00215           aMinuit->mnparm(0, "X0 (cm)",pval[0],0.01,-2,2,error_flag);
00216           aMinuit->mnparm(1, "Y0 (cm)",pval[1],0.01,-2,2,error_flag);
00217           aMinuit->mnparm(2, "Ux",pval[2],0.001,-0.02,0.02,error_flag);
00218           aMinuit->mnparm(3, "Uy",pval[3],0.001,-0.02,0.02,error_flag);
00219           arglist[0] = 500.;
00220           arglist[1] = 1.;
00221           if (mode < 3)
00222             aMinuit->mnexcm(com,arglist,2,error_flag);  
00223           else
00224             {
00225               aMinuit->mnexcm(com1,arglist,2,error_flag);
00226               aMinuit->mnexcm(com2,arglist,2,error_flag);
00227               aMinuit->mnexcm(com1,arglist,2,error_flag);
00228               aMinuit->mnexcm(com2,arglist,2,error_flag);
00229             }
00230           aMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
00231           aMinuit->mnprin(3,amin);
00232           for(int ip=0;ip<mxPar;ip++) {
00233             aMinuit->mnpout(ip,paraName[ip],pval[ip],perr[ip],plo[ip],phi[ip],istat);}
00234           printf("PARAValues X,Y,UX,UY for %s = %f , %f , %f, %f  \n",core,pval[0],pval[1],pval[2],pval[3]);
00235           printf("PARAErrors X,Y,UX,UY for %s = %f , %f , %f, %f  \n",core,perr[0],perr[1],perr[2],perr[3]);
00236           aMinuit->Delete();
00237         }
00238       }
00239       
00240     TVector3 V(pval[0],pval[1],0);
00241     TVector3 U(pval[2],pval[3],1);
00242     float errxy = 0;
00243     float errz = 0;
00244     int nTrackxy=0;
00245     int nTrackz=0;
00246     float chi2 = 0;
00247     U=U.Unit();
00248     if (1){
00249       vector<TrackStump>::iterator it;
00250       for (it=util.track.begin() ; it < util.track.end(); it++) {
00251         TrackStump *t= &(*it);  
00252         if (t->bad) continue;
00253         double lambda=0;double kappa=0;
00254         TVector3 r(t->r),p(t->p);
00255         twoLineDca3D(lambda,kappa,V,U,r,p);
00256         TVector3 D=((r-V)+(lambda*p-kappa*U));
00257         double dT2=D.Perp2(); // compute distance only once
00258         double dZ2=D.z()*D.z();
00259         if((dT2+t->ery2)<dmax2){
00260           nTrackxy++;
00261           errxy+=t->ery2;
00262           // chi2+=dT2/t->ery2;
00263           chi2+=exp(-1*D.x()*D.x()/(2*t->ery2))+exp(-1*D.y()*D.y()/(2*t->ery2));
00264         }
00265         if((dZ2+t->erz2)<dmax2){
00266           errz+=t->erz2;nTrackz++;
00267           // chi2+=dZ2/t->erz2;
00268           chi2+=exp(-1*dZ2/(2*t->erz2));
00269         }
00270         //if((dT2+t->ery2<dmax2)&&(dZ2+t->erz2<dmax2))
00271         // chi2+=exp(-1*D.Mag()*D.Mag()/(2*t->ery2+2*t->erz2));
00272       }
00273       
00274     }
00275     printf("Uncertainty for %s = %i , %f ,%i , %f chi^2 %f  \n",core,nTrackxy,errxy,nTrackz,errz,chi2);
00276     gMinuit->Delete();
00277     }
00278     util.hA[19]->GetXaxis()->SetTitle(core);
00279     util.hA[19]->GetYaxis()->SetTitle("cm");
00280     
00281     printf("\n    ***done!***, scan Chi2 around solution\n");    
00282     util.evalSolution(pval);
00283     util.scanChi2(pval,0); //0=Vx-Vy
00284 #ifdef FIT_XYSLOPES
00285     util.scanChi2(pval,1); //1=Vx-Ux
00286 #endif  
00287     // Save histograms
00288     HList.ls();
00289     // return 0;
00290     TString outName="out/3D_beam_"; outName+=core;outName+=".hist.root";
00291     TFile f( outName,"recreate");
00292     assert(f.IsOpen());
00293     printf("%d histos are written  to '%s' ...\n",HList.GetEntries(),outName.Data());
00294     HList.Write();
00295     f.Close();
00296     //   gMinuit->Delete();
00297 
00298   }
00299 }
00300 
00301