StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFgtClustEvalMaker.cxx
1 // *-- Author : J.Balewski
2 //
3 // $Id: StFgtClustEvalMaker.cxx,v 1.2 2011/04/11 19:35:38 fisyak Exp $
4 
5 #include <TVector3.h>
6 #include <TH2.h>
7 #include <TF1.h>
8 #include <TFile.h>
9 #include <TLine.h>
10 #include <TPolyLine.h>
11 #include <TCrown.h>
12 #include <TRandom3.h>
13 #include "TMath.h"
14 #include "StFgtClustEvalMaker.h"
15 #include "StFgtSlowSimuMaker.h"
16 #include "StFgtClustFindMaker.h"
17 
18 ClassImp(StFgtClustEvalMaker)
19 
20 //--------------------------------------------
21 StFgtClustEvalMaker::StFgtClustEvalMaker(const char *name):StMaker(name){
23  setHList(0);
24  memset(hA,0,sizeof(hA));
25  geom=new StFgtGeom();
26  par_minDelRad=0.07; //(cm) radial match cut off ,default=0.07
27  par_minRdPhi=0.02; //(cm) match cut off ,default=0.02
28 }
29 
30 //--------------------------------------------
31 void
32 StFgtClustEvalMaker::Clear(Option_t *) {
33  StMaker::Clear();
34 }
35 
36 //--------------------------------------------
37 StFgtClustEvalMaker::~StFgtClustEvalMaker(){
38 
39 }
40 
41 //_______________________________________________
42 //________________________________________________
43 void
44 StFgtClustEvalMaker::saveHisto(TString fname){
45  TString outName=fname+".hist.root";
46  TFile f( outName,"recreate");
47  assert(f.IsOpen());
48  printf("%d histos are written to '%s' ...\n",HList->GetEntries(),outName.Data());
49 
50  HList->Write();
51  f.Close();
52 
53 }
54 
55 
56 //--------------------------------------------
57 //--------------------------------------------
58 //--------------------------------------------
59 Int_t
60 StFgtClustEvalMaker::Finish(){
61  LOG_INFO<<"::Finish() \n"<< endm;
62 
63 
64  return StMaker::Finish();
65 }
66 
67 //--------------------------------------------
68 //--------------------------------------------
69 //--------------------------------------------
70 Int_t
71 StFgtClustEvalMaker::Make(){
72  mInpEve++;
73 
74  LOG_INFO <<"::Make() inpEve="<< mInpEve<<endm;
75  StFgtSlowSimuMaker *ssMk=(StFgtSlowSimuMaker *) GetMaker("FgtSlowSimu");
76  assert(ssMk);
77  StFgtClustFindMaker *cfMk=(StFgtClustFindMaker *) GetMaker("FgtClustFind");
78  assert(cfMk);
79 
80  printf(" INPUT: \n disk #g2t_tracks #Rad_clust #Phi_clust \n");
81 
82 
83  // .... match 1D clusters in every plane of all disks
84  for(int iDisk=0;iDisk<kFgtMxDisk;iDisk++) {
85  int nG=0;
86  for(int iQuad=0; iQuad<kFgtMxQuad;iQuad++) nG+=ssMk->mG2tHitList[iDisk][iQuad].size();
87 
88  int nrR= cfMk->mRadClustList[iDisk].size();
89  int nrP= cfMk->mPhiClustList[iDisk].size();
90 
91  int nmR=matchRadClust1D( ssMk->mG2tHitList[iDisk],cfMk->mRadClustList[iDisk]);
92  int nmP=matchPhiClust1D( ssMk->mG2tHitList[iDisk],cfMk->mPhiClustList[iDisk]);
93  // printf("match Disk=%d nG=%d nrR=%d->%d nrP=%d->%d\n",iDisk+1,nG,nrR,nmR, nrP,nmP);
94 
95  hA[0]->Fill(10*iDisk+0,nG);
96  hA[0]->Fill(10*iDisk+1,nrR);
97  hA[0]->Fill(10*iDisk+2,nrP);
98  hA[0]->Fill(10*iDisk+3,nmR);
99  hA[0]->Fill(10*iDisk+4,nmP);
100  if(nrR>nG) printf("#OV Rad %d %d\n",nrR,nG);
101  if(nrR>nmR) printf("#MI Rad %d %d\n",nrR,nmR);
102  if(nrP>nG) printf("#OV Phi %d %d\n",nrP,nG);
103  if(nrP>nmP) printf("#MI Phi %d %d\n",nrP,nmP);
104  }
105  return kStOK;
106 }
107 
108 
109 //--------------------------------------------
110 //--------------------------------------------
111 //--------------------------------------------
112 int
113 StFgtClustEvalMaker::matchRadClust1D( vector<fgt_g2t_auxil> *g2tTrL, vector<fgt_cluster1D> &clL){
114 
115  /* brute force used,
116  unnecessary scans over all clusters in one disk, could do it by quadrant
117  */
118 
119  int nM=0;
120  for(int iQuad=0; iQuad<kFgtMxQuad;iQuad++)
121  for(UInt_t ih=0;ih<g2tTrL[iQuad].size();ih++) {
122  double gR=g2tTrL[iQuad][ih].Rlab.Perp(); // radial distance from Z-axis
123  double minDelR=par_minDelRad;
124  int jc=-1;// best match
125  for(UInt_t ic=0;ic<clL.size();ic++) {
126  if(clL[ic].matched) continue; //match only once
127  double delR=clL[ic].position - gR;
128  if(TMath::Abs(minDelR)<TMath::Abs(delR)) continue;
129  minDelR=delR;// better match found
130  jc=ic;
131  }
132  if(jc>=0) {
133  hA[5]->Fill(minDelR*1e4); // in um
134  hA[7]->Fill(gR,minDelR*1e4); // in um
135  hA[9]->Fill(g2tTrL[iQuad][ih].Rlab.x(),g2tTrL[iQuad][ih].Rlab.y());
136  // printf("match Rad: iQ=%d gT=%d gR=%.2f match jc=%d minDelR=%.3f\n",iQuad, ih,gR,jc,minDelR);
137  nM++;
138  clL[jc].matched=true;
139  } //else printf("nomatch ih=%d gR=%.2f\n", ih,gR);
140 
141  }
142  return nM;
143 }
144 
145 //--------------------------------------------
146 //--------------------------------------------
147 //--------------------------------------------
148 int
149 StFgtClustEvalMaker::matchPhiClust1D( vector<fgt_g2t_auxil> *g2tTrL, vector<fgt_cluster1D> &clL){
150  const double pi=2.*acos(0.), dpi=2.*pi; // double pi
151 
152  /* brute force used,
153  unnecessary scans over all clusters in one disk, could do it by quadrant
154  */
155 
156  int nM=0;
157  for(int iQuad=0; iQuad<kFgtMxQuad;iQuad++)
158  for(UInt_t ih=0;ih<g2tTrL[iQuad].size();ih++) {
159  double gR=g2tTrL[iQuad][ih].Rlab.Perp(); // radial distance from Z-axis
160  double gPhi=g2tTrL[iQuad][ih].Rlab.Phi(); // radial distance from Z-axis
161  while (gPhi<-pi) gPhi+=dpi;
162  while (gPhi>pi) gPhi-=dpi; // angle in [0,2pi]
163 
164  double minRdPhi=par_minRdPhi;
165  int jc=-1;// best match
166  for(UInt_t ic=0;ic<clL.size();ic++) {
167  if(clL[ic].matched) continue; //match only once
168  double rPhi=clL[ic].position;
169  while (rPhi<-pi) rPhi+=dpi;
170  while (rPhi>pi) rPhi-=dpi;// angle in [0,2pi]
171  double delPhi=gPhi-rPhi;
172  if(delPhi>pi) delPhi-=dpi;
173  if(delPhi<-pi) delPhi+=dpi;
174  double RdelPhi=gR*delPhi;
175  //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);
176 
177  if(TMath::Abs(minRdPhi)<TMath::Abs(RdelPhi)) continue;
178  minRdPhi=RdelPhi;// better match found
179  jc=ic;
180  }
181  if(jc>=0) {
182  hA[6]->Fill(minRdPhi*1e4); // in um
183  hA[8]->Fill(gR,minRdPhi*1e4); // in um
184  hA[10]->Fill(g2tTrL[iQuad][ih].Rlab.x(),g2tTrL[iQuad][ih].Rlab.y());
185  //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);
186  nM++;
187  clL[jc].matched=true;
188  } //else printf("nomatch ih=%d gR=%.2f gPhi/deg=%.2f\n", ih,gR,gPhi/3.1417*180.);
189 
190  }
191  return nM;
192 }
193 
194 //--------------------------------------------
195 //--------------------------------------------
196 //--------------------------------------------
197 Int_t
198 StFgtClustEvalMaker::Init(){
199  assert(HList);
200  mInpEve=0;
201 // int i;
202 
203  hA[0]=new TH1F("ev_Stat1D","Eval clust matching; x=10*disk",100, -0.5,99.5);
204 
205  double Yrange=2*geom->radStrip_pitch()*1e4;
206 
207  hA[5]= new TH1F("ev_errRad","error in reco R (track-recoClust); R_err (um)",100,-par_minDelRad*1e4,par_minDelRad*1e4);
208  hA[6]= new TH1F("ev_RerrPhi","error in reco phi*R (track-recoClust); R*#Phi_err (um)",100,-par_minRdPhi*1e4,par_minRdPhi*1e4);
209 
210 
211 
212  hA[7]=new TH2F("ev_dRad_R","#DeltaR accuracy of reco cluster ;thrown R (cm); #DeltaR (um)",45,0.,45.,50,-Yrange,Yrange);
213  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);
214 
215 
216  char text[1000];
217  sprintf(text,"gen clust X-Y , R-matched (#DeltaR<%.3fmm);Lab X (cm); Lab Y (cm); ",par_minDelRad*10);
218  hA[9]=new TH2F("ev_RokXY",text,90,-45.,45.,90,-45.,45.);
219 
220  sprintf(text,"gen clust X-Y , #Phi-matched (R#Delta#Phi<%.3fmm);Lab X (cm); Lab Y (cm); ",par_minRdPhi*10);
221  hA[10]=new TH2F("ev_PhiokXY",text,90,-45.,45.,90,-45.,45.);
222 
223 
224  //..............add gadgets to histos.....................
225  TList *Lx; TLine *ln; TPolyLine *pln;
226 
227  //... strip size for cluster error plots: delRad
228  double x[100],y[100];
229 
230  Lx=hA[7]->GetListOfFunctions(); // delR
231  double y1=par_minDelRad*1e4;
232  ln=new TLine(0,y1,45,y1); ln->SetLineColor(kGreen);Lx->Add(ln);
233  ln=new TLine(0,-y1,45,-y1); ln->SetLineColor(kGreen);Lx->Add(ln);
234  x[0]=x[9]=geom->Rin(); y[0]=y[9]=0;
235  x[1]=x[8]=x[0]; y[1]=geom->radStrip_pitch()*1e4/2; y[8]=-y[1];
236  x[2]=x[7]=geom->Rmid(); y[2]=y[1]; y[7]=-y[2];
237  x[3]=x[6]=x[2]; y[3]=geom->radStrip_pitch()*1e4/2; y[6]=-y[3];
238  x[4]=x[5]=geom->Rout(); y[4]=y[3]; y[5]=-y[4];
239  pln=new TPolyLine(10,x,y);
240  pln->SetLineColor(kRed); Lx->Add(pln);
241 
242 
243  Lx=hA[8]->GetListOfFunctions(); // R*delPhi
244  y1=par_minRdPhi*1e4;
245  ln=new TLine(0,y1,45,y1); ln->SetLineColor(kGreen);Lx->Add(ln);
246  ln=new TLine(0,-y1,45,-y1); ln->SetLineColor(kGreen);Lx->Add(ln);
247  double hl=geom->phiStrip_pitch()*1e4/2 * geom->Rin();
248  double hh=2*hl;
249  double hm=hh*geom->Rmid()/geom->Rmid();
250  x[0]=x[12]=x[11]=geom->Rin(); y[0]=y[12]=hl; y[11]=-y[0];
251  x[1]=x[10]=geom->Rmid(); y[1]=hh; y[10]=-y[1];
252  x[2]=x[9]=x[1]; y[2]=hm; y[9]=-y[2];
253  x[3]=x[8]=geom->Rmid(); y[3]=hh; y[8]=-y[3];
254  x[4]=x[7]=x[3]; y[4]=hl; y[7]=-y[4];
255  x[5]=x[6]=geom->Rout(); y[5]=hh; y[6]=-y[5];
256  pln=new TPolyLine(13,x,y); pln->SetLineColor(kRed);Lx->Add(pln);
257 
258 
259  for(int i=0;i<mxH;i++)
260  if(hA[i]) HList->Add(hA[i]);
261 
262  LOG_INFO<<Form("::Init params match limit (cm) delR=%.3f R*delPhi=%.3f", par_minDelRad,par_minRdPhi)<< endm;
263 
264 
265  return StMaker::Init();
266 }
267 
268 
271 
272 // $Log: StFgtClustEvalMaker.cxx,v $
273 // Revision 1.2 2011/04/11 19:35:38 fisyak
274 // Replace uint by UInt_t, use TMath
275 //
276 // Revision 1.1 2011/04/07 19:31:22 balewski
277 // start
278 //
279 
280 
281 
282 
283 
StFgtSlowSimuMaker::mG2tHitList
vector< fgt_g2t_auxil > mG2tHitList[kFgtMxDisk+1][kFgtMxQuad]
accepted track segements
Definition: StFgtSlowSimuMaker.h:75
StMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
kStOK
Definition: Stypes.h:40
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
StFgtClustEvalMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StFgtClustEvalMaker.cxx:32
Generated by   doxygen 1.8.5