StRoot: StEEmcPool/FgtSandbox1/StFgtClustEvalMaker.cxx Source File

00001 // *-- Author : J.Balewski
00002 // 
00003 // $Id: StFgtClustEvalMaker.cxx,v 1.2 2011/04/11 19:35:38 fisyak Exp $
00004 
00005 #include <TVector3.h>
00006 #include <TH2.h>
00007 #include <TF1.h>
00008 #include <TFile.h>
00009 #include <TLine.h>
00010 #include <TPolyLine.h>
00011 #include <TCrown.h>
00012 #include <TRandom3.h>
00013 #include "TMath.h"
00014 #include "StFgtClustEvalMaker.h"
00015 #include "StFgtSlowSimuMaker.h"
00016 #include "StFgtClustFindMaker.h"
00017 
00018 ClassImp(StFgtClustEvalMaker)
00019 
00020 //--------------------------------------------
00021 StFgtClustEvalMaker::StFgtClustEvalMaker(const char *name):StMaker(name){
00023   setHList(0);
00024   memset(hA,0,sizeof(hA));
00025   geom=new StFgtGeom();
00026   par_minDelRad=0.07; //(cm) radial match cut off ,default=0.07
00027   par_minRdPhi=0.02; //(cm)  match cut off ,default=0.02
00028 }
00029 
00030 //--------------------------------------------
00031 void 
00032 StFgtClustEvalMaker::Clear(Option_t *) {
00033    StMaker::Clear();
00034 }
00035 
00036 //--------------------------------------------
00037 StFgtClustEvalMaker::~StFgtClustEvalMaker(){
00038 
00039 }
00040 
00041 //_______________________________________________
00042 //________________________________________________
00043 void
00044 StFgtClustEvalMaker::saveHisto(TString fname){
00045   TString outName=fname+".hist.root";
00046   TFile f( outName,"recreate");
00047   assert(f.IsOpen());
00048   printf("%d histos are written  to '%s' ...\n",HList->GetEntries(),outName.Data());
00049 
00050   HList->Write();
00051   f.Close();
00052 
00053 }
00054 
00055 
00056 //--------------------------------------------
00057 //--------------------------------------------
00058 //--------------------------------------------
00059 Int_t 
00060 StFgtClustEvalMaker::Finish(){
00061   LOG_INFO<<"::Finish() \n"<<  endm; 
00062 
00063 
00064   return StMaker::Finish();
00065 }
00066 
00067 //--------------------------------------------
00068 //--------------------------------------------
00069 //--------------------------------------------
00070 Int_t 
00071 StFgtClustEvalMaker::Make(){
00072   mInpEve++;
00073   
00074   LOG_INFO <<"::Make() inpEve="<<  mInpEve<<endm;
00075   StFgtSlowSimuMaker *ssMk=(StFgtSlowSimuMaker *) GetMaker("FgtSlowSimu");
00076   assert(ssMk);
00077   StFgtClustFindMaker *cfMk=(StFgtClustFindMaker *) GetMaker("FgtClustFind");
00078   assert(cfMk);
00079   
00080   printf(" INPUT: \n disk  #g2t_tracks  #Rad_clust   #Phi_clust \n");
00081 
00082 
00083   // .... match 1D clusters in every plane of all disks
00084   for(int iDisk=0;iDisk<kFgtMxDisk;iDisk++) {
00085     int nG=0;
00086     for(int iQuad=0; iQuad<kFgtMxQuad;iQuad++) nG+=ssMk->mG2tHitList[iDisk][iQuad].size();
00087 
00088     int nrR= cfMk->mRadClustList[iDisk].size();
00089     int nrP= cfMk->mPhiClustList[iDisk].size();
00090 
00091     int nmR=matchRadClust1D( ssMk->mG2tHitList[iDisk],cfMk->mRadClustList[iDisk]);
00092     int nmP=matchPhiClust1D( ssMk->mG2tHitList[iDisk],cfMk->mPhiClustList[iDisk]);
00093     // printf("match Disk=%d  nG=%d   nrR=%d->%d  nrP=%d->%d\n",iDisk+1,nG,nrR,nmR, nrP,nmP);
00094 
00095     hA[0]->Fill(10*iDisk+0,nG);
00096     hA[0]->Fill(10*iDisk+1,nrR);
00097     hA[0]->Fill(10*iDisk+2,nrP);
00098     hA[0]->Fill(10*iDisk+3,nmR);
00099     hA[0]->Fill(10*iDisk+4,nmP);
00100     if(nrR>nG) printf("#OV Rad %d %d\n",nrR,nG);
00101     if(nrR>nmR) printf("#MI Rad %d %d\n",nrR,nmR);
00102     if(nrP>nG) printf("#OV Phi %d %d\n",nrP,nG);
00103     if(nrP>nmP) printf("#MI Phi %d %d\n",nrP,nmP);
00104   }
00105  return kStOK;
00106 }
00107 
00108  
00109 //--------------------------------------------
00110 //--------------------------------------------
00111 //--------------------------------------------
00112 int
00113 StFgtClustEvalMaker::matchRadClust1D( vector<fgt_g2t_auxil> *g2tTrL,  vector<fgt_cluster1D> &clL){
00114 
00115   /* brute force used,
00116      unnecessary scans over all clusters in one disk, could do it by quadrant
00117   */
00118 
00119   int nM=0;
00120   for(int iQuad=0; iQuad<kFgtMxQuad;iQuad++)     
00121     for(UInt_t ih=0;ih<g2tTrL[iQuad].size();ih++) {
00122       double gR=g2tTrL[iQuad][ih].Rlab.Perp(); // radial distance from Z-axis
00123       double minDelR=par_minDelRad;
00124       int jc=-1;// best match
00125       for(UInt_t ic=0;ic<clL.size();ic++) {
00126         if(clL[ic].matched) continue; //match only once
00127         double delR=clL[ic].position - gR;
00128         if(TMath::Abs(minDelR)<TMath::Abs(delR)) continue;
00129         minDelR=delR;// better match found
00130         jc=ic;
00131       }
00132       if(jc>=0) { 
00133         hA[5]->Fill(minDelR*1e4); // in um
00134         hA[7]->Fill(gR,minDelR*1e4); // in um
00135         hA[9]->Fill(g2tTrL[iQuad][ih].Rlab.x(),g2tTrL[iQuad][ih].Rlab.y());
00136         //      printf("match Rad: iQ=%d gT=%d gR=%.2f  match  jc=%d  minDelR=%.3f\n",iQuad, ih,gR,jc,minDelR);
00137         nM++;
00138         clL[jc].matched=true;
00139       } //else  printf("nomatch ih=%d gR=%.2f\n", ih,gR);
00140 
00141     }
00142   return nM;
00143 }
00144 
00145 //--------------------------------------------
00146 //--------------------------------------------
00147 //--------------------------------------------
00148 int
00149 StFgtClustEvalMaker::matchPhiClust1D( vector<fgt_g2t_auxil> *g2tTrL,  vector<fgt_cluster1D> &clL){
00150   const double pi=2.*acos(0.),  dpi=2.*pi; // double pi
00151 
00152   /* brute force used,
00153      unnecessary scans over all clusters in one disk, could do it by quadrant
00154   */
00155 
00156   int nM=0;
00157   for(int iQuad=0; iQuad<kFgtMxQuad;iQuad++)     
00158     for(UInt_t ih=0;ih<g2tTrL[iQuad].size();ih++) {
00159       double gR=g2tTrL[iQuad][ih].Rlab.Perp(); // radial distance from Z-axis
00160       double gPhi=g2tTrL[iQuad][ih].Rlab.Phi(); // radial distance from Z-axis
00161       while (gPhi<-pi) gPhi+=dpi;
00162       while (gPhi>pi) gPhi-=dpi; // angle in [0,2pi]
00163 
00164       double minRdPhi=par_minRdPhi;
00165       int jc=-1;// best match
00166       for(UInt_t ic=0;ic<clL.size();ic++) {
00167         if(clL[ic].matched) continue; //match only once
00168         double rPhi=clL[ic].position;
00169         while (rPhi<-pi) rPhi+=dpi;
00170         while (rPhi>pi) rPhi-=dpi;// angle in [0,2pi]
00171         double delPhi=gPhi-rPhi;
00172         if(delPhi>pi) delPhi-=dpi;
00173         if(delPhi<-pi) delPhi+=dpi;
00174         double RdelPhi=gR*delPhi;
00175         //printf("phi1,2: g r=%.3f phi/deg=%.2f, r phi/deg=%.2f,  delPR/cm=%.3f  \n",gR,gPhi/3.1416*180.,rPhi/3.1416*180.,RdelPhi);
00176 
00177         if(TMath::Abs(minRdPhi)<TMath::Abs(RdelPhi)) continue;
00178         minRdPhi=RdelPhi;// better match found
00179         jc=ic;
00180       }
00181       if(jc>=0) { 
00182         hA[6]->Fill(minRdPhi*1e4); // in um
00183         hA[8]->Fill(gR,minRdPhi*1e4); // in um
00184         hA[10]->Fill(g2tTrL[iQuad][ih].Rlab.x(),g2tTrL[iQuad][ih].Rlab.y());
00185         //printf("match Phi iQ=%d gT=%d gR=%.2f gPhi/deg=%.2f  match  jc=%d  minRdPhi=%.3f\n",iQuad, ih,gR,gPhi/3.1417*180.,jc,minRdPhi);
00186         nM++;
00187         clL[jc].matched=true;
00188       } //else  printf("nomatch ih=%d gR=%.2f gPhi/deg=%.2f\n", ih,gR,gPhi/3.1417*180.);
00189 
00190     }
00191   return nM;
00192 }
00193 
00194 //--------------------------------------------
00195 //--------------------------------------------
00196 //--------------------------------------------
00197 Int_t 
00198 StFgtClustEvalMaker::Init(){
00199   assert(HList);
00200   mInpEve=0;
00201 //  int i;
00202 
00203   hA[0]=new TH1F("ev_Stat1D","Eval clust matching; x=10*disk",100, -0.5,99.5);
00204 
00205   double Yrange=2*geom->radStrip_pitch()*1e4;
00206   
00207   hA[5]= new TH1F("ev_errRad","error in reco  R (track-recoClust); R_err (um)",100,-par_minDelRad*1e4,par_minDelRad*1e4);
00208   hA[6]= new TH1F("ev_RerrPhi","error in reco phi*R   (track-recoClust); R*#Phi_err (um)",100,-par_minRdPhi*1e4,par_minRdPhi*1e4);
00209 
00210 
00211 
00212   hA[7]=new TH2F("ev_dRad_R","#DeltaR    accuracy of reco cluster ;thrown R (cm); #DeltaR   (um)",45,0.,45.,50,-Yrange,Yrange);
00213   hA[8]=new TH2F("ev_RdPhi_R","R*#Delta#Phi   accuracy of reco cluster ;thrown R (cm);R*#Delta#Phi (um)",45,0.,45.,200,-Yrange,Yrange);
00214 
00215 
00216   char text[1000];
00217   sprintf(text,"gen clust X-Y , R-matched (#DeltaR<%.3fmm);Lab  X (cm); Lab Y (cm); ",par_minDelRad*10);
00218   hA[9]=new TH2F("ev_RokXY",text,90,-45.,45.,90,-45.,45.);
00219 
00220   sprintf(text,"gen clust X-Y , #Phi-matched (R#Delta#Phi<%.3fmm);Lab  X (cm); Lab Y (cm); ",par_minRdPhi*10);
00221   hA[10]=new TH2F("ev_PhiokXY",text,90,-45.,45.,90,-45.,45.);
00222 
00223 
00224   //..............add gadgets to histos.....................
00225   TList *Lx; TLine *ln;    TPolyLine *pln;  
00226 
00227   //... strip size for cluster error plots: delRad
00228   double x[100],y[100];
00229 
00230    Lx=hA[7]->GetListOfFunctions(); // delR
00231    double y1=par_minDelRad*1e4;
00232    ln=new TLine(0,y1,45,y1); ln->SetLineColor(kGreen);Lx->Add(ln);
00233    ln=new TLine(0,-y1,45,-y1); ln->SetLineColor(kGreen);Lx->Add(ln);
00234    x[0]=x[9]=geom->Rin();  y[0]=y[9]=0;
00235    x[1]=x[8]=x[0];         y[1]=geom->radStrip_pitch()*1e4/2; y[8]=-y[1];
00236    x[2]=x[7]=geom->Rmid(); y[2]=y[1];                y[7]=-y[2];
00237    x[3]=x[6]=x[2];         y[3]=geom->radStrip_pitch()*1e4/2; y[6]=-y[3];
00238    x[4]=x[5]=geom->Rout(); y[4]=y[3];                y[5]=-y[4];
00239    pln=new TPolyLine(10,x,y);
00240    pln->SetLineColor(kRed); Lx->Add(pln);
00241 
00242 
00243    Lx=hA[8]->GetListOfFunctions(); // R*delPhi
00244    y1=par_minRdPhi*1e4;
00245    ln=new TLine(0,y1,45,y1); ln->SetLineColor(kGreen);Lx->Add(ln);
00246    ln=new TLine(0,-y1,45,-y1); ln->SetLineColor(kGreen);Lx->Add(ln);
00247    double hl=geom->phiStrip_pitch()*1e4/2 * geom->Rin();
00248    double hh=2*hl;
00249    double hm=hh*geom->Rmid()/geom->Rmid();
00250    x[0]=x[12]=x[11]=geom->Rin();  y[0]=y[12]=hl;  y[11]=-y[0];
00251    x[1]=x[10]=geom->Rmid();  y[1]=hh;     y[10]=-y[1];
00252    x[2]=x[9]=x[1];           y[2]=hm;     y[9]=-y[2];
00253    x[3]=x[8]=geom->Rmid(); y[3]=hh;     y[8]=-y[3];
00254    x[4]=x[7]=x[3];           y[4]=hl;     y[7]=-y[4];
00255    x[5]=x[6]=geom->Rout();   y[5]=hh;     y[6]=-y[5];
00256    pln=new TPolyLine(13,x,y);  pln->SetLineColor(kRed);Lx->Add(pln);
00257 
00258 
00259   for(int i=0;i<mxH;i++) 
00260     if(hA[i]) HList->Add(hA[i]);
00261   
00262   LOG_INFO<<Form("::Init params match limit (cm) delR=%.3f  R*delPhi=%.3f", par_minDelRad,par_minRdPhi)<<  endm; 
00263 
00264 
00265   return StMaker::Init();
00266 }
00267 
00268 
00271 
00272 // $Log: StFgtClustEvalMaker.cxx,v $
00273 // Revision 1.2  2011/04/11 19:35:38  fisyak
00274 // Replace uint by UInt_t, use TMath
00275 //
00276 // Revision 1.1  2011/04/07 19:31:22  balewski
00277 // start
00278 //
00279 
00280 
00281  
00282 
00283