StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFtpcLaserCalib.cc
1 // $Id: StFtpcLaserCalib.cc,v 1.11 2010/01/13 15:43:33 jcs Exp $
2 //
3 // $Log: StFtpcLaserCalib.cc,v $
4 // Revision 1.11 2010/01/13 15:43:33 jcs
5 // improved LOG_DEBUG output
6 //
7 // Revision 1.10 2009/12/09 14:41:30 jcs
8 // delta_t0 and delta_gas can now both = 0
9 // new space point calculation always necessary since reconstruction done with data t0
10 //
11 // Revision 1.9 2009/10/14 15:59:55 jcs
12 // changes to be able to vary the gas temperature in addition to varying t0 and
13 // gas composition
14 //
15 // Revision 1.8 2009/10/06 14:51:45 jcs
16 // exit laser_fit to avoid FPE if either <=2 hits on track or if helix fit fails to converge
17 //
18 // Revision 1.7 2008/05/15 22:39:47 jcs
19 // re-activate helix fit
20 //
21 // Revision 1.6 2007/01/22 13:08:15 jcs
22 // replace cout with LOG
23 //
24 // Revision 1.5 2007/01/19 08:22:59 jcs
25 // replace true, false with kTRUE, kFALSE
26 //
27 // Revision 1.4 2006/06/19 12:37:58 jcs
28 // Use negative value of zfieldcage for FTPC east
29 //
30 // Revision 1.3 2006/04/04 10:57:05 jcs
31 // Fix memory leak
32 //
33 // Revision 1.2 2006/03/15 15:13:56 jcs
34 // add lines for listing CVS update info
35 //
36 
37 #include "StFtpcLaserCalib.hh"
38 #include "StFtpcTrackMaker/StFtpcPoint.hh"
39 #include "StFtpcTrackMaker/StFtpcTrack.hh"
40 #include "StFtpcTrackMaker/StFtpcVertex.hh"
41 #include "TObjArray.h"
42 
43 #define HELIX_FIT
44 
45 const Float_t rad2grad=180/TMath::Pi();
46 
47 // define functions etc. used for fit added on 11/21/2003
48 static Double_t funcxz(float x,float z,Double_t *par);
49 static Double_t funcyz(float y,float z,Double_t *par);
50 static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
51 
52 //
53 
54 // im prinzip unschoen !!! Warning : defined but not used !!!
55 static int nhits; // needed by the fcn used for minuit.
56 static Float_t z[11],x[11],y[11];
57 static Float_t ex[11],ey[11];
58 
59 //______________________________________________________________________________
60 
61 StFtpcLaserCalib::StFtpcLaserCalib()
62 {
63  mMinuit=new TMinuit(4);
64  //mMinuit->SetPrintLevel(-1);
65 }
66 
67 //______________________________________________________________________________
68 
69 //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)
70 StFtpcLaserCalib::StFtpcLaserCalib(int ftpc, int lsec, int straight, int gfit, int minz, int maxz, int minrad, int maxrad, float gt0, float ggas, float gTemp, StFtpcLaserTrafo *trafo, StarMagField *gmagf)
71 {
72  FTPC=ftpc;
73  LSEC=lsec;
74  STRAIGHT=straight;
75  MINZ=minz;
76  MAXZ=maxz;
77  MINRAD=minrad;
78  MAXRAD=maxrad;
79  GAUSFIT=gfit;
80  mMinuit=new TMinuit(4);
81  deltat0=gt0;
82  deltagas=ggas;
83  deltaTemp=gTemp;
84 
85  usedfit=0;
86 
87 #ifdef HELIX_FIT
88  trcharge=0;
89  p=pt=invp=invpt=0;
90 #endif
91 
92  magf=gmagf;
93 
94  //if (deltat0!=0 || deltagas!=0)
95  mtrafo=trafo;
96 }
97 
98 //______________________________________________________________________________
99 
100 //static
101 Double_t funcxz(float x,float z,Double_t *par)
102 {
103 
104  return x-(par[0]*z+par[1]);
105 }
106 
107 //static
108 Double_t funcyz(float y,float z,Double_t *par)
109 {
110  return y-(par[2]*z+par[3]);
111 }
112 
113 //______________________________________________________________________________
114 
115 //static
116 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
117 //void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
118 {
119 
120  Int_t i; //calculate chisquare
121 
122  Double_t chisq = 0;
123  Double_t deltaxz, deltayz;
124  //LOG_DEBUG<<"in fcn vor for .."<<nhits<<" hits"<<endm;
125 
126  for (i=0;i<nhits; i++)
127  {
128  // Fehler aus Daten noch intgerieren !
129  // so richtig !??????
130  deltaxz=funcxz(x[i],z[i],par)/ex[i];
131  deltayz=funcyz(y[i],z[i],par)/ey[i];
132  chisq += (deltaxz*deltaxz+deltayz*deltayz);
133  }
134 
135  //LOG_DEBUG<<" f = "<<f<<endm;
136  f = chisq;
137 }
138 
139 //______________________________________________________________________________
140 
141 bool StFtpcLaserCalib::cut(int i)
142 {
143  LOG_DEBUG<<"StFtpcLaserCalib::cut i "<<i<<" z "<< z[i]<<" MINZ "<<MINZ<<" MAXZ "<<MAXZ<<" radius "<<radius[i]<<" MINRAD "<<MINRAD<<" MAXRAD "<<MAXRAD<< " laser_sector "<<hsec[i]<<endm;
144  if ((z[i]>MINZ && z[i]<MAXZ) && (radius[i]>MINRAD && radius[i]<MAXRAD) && laser_sector(FTPC,LSEC,hsec[i]) )
145  return kTRUE;
146  else
147  return kFALSE;
148 }
149 
150 //______________________________________________________________________________
151 
152 void StFtpcLaserCalib::minuit_init()
153 {
154  //LOG_DEBUG<<"Minuit init ..."<<endm;
155  //mMinuit = new TMinuit(4); //???? why is there a link error !?
156  mMinuit->SetPrintLevel(-1);
157  mMinuit->SetFCN(&fcn);
158  //LOG_DEBUG<<"set fcn"<<endm;
159  ierflg = 0;
160  arglist[0] = 1;
161 
162  mMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
163  //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm;
164 }
165 
166 //______________________________________________________________________________
167 
168 void StFtpcLaserCalib::minuit_set_par()
169 {
170  mMinuit->SetPrintLevel(-1);
171  //LOG_DEBUG<<"Minuit set parameter ..."<<endm;
172  //static Double_t vstart[4] = {1,1,1,1};
173  Double_t vstart[4] = {1,1,1,1};
174  // ggf. daten werte als Startparameter benutzen !
175  //static Double_t step[4] = {0.1,0.1,0.1,0.1};
176  Double_t step[4] = {0.1,0.1,0.1,0.1};
177  mMinuit->mnparm(0, "mxz", vstart[0], step[0], 0,0,ierflg);
178  mMinuit->mnparm(1, "bxz", vstart[1], step[1], 0,0,ierflg);
179  mMinuit->mnparm(2, "myz", vstart[2], step[2], 0,0,ierflg);
180  mMinuit->mnparm(3, "byz", vstart[3], step[3], 0,0,ierflg);
181  //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm;
182 }
183 
184 //______________________________________________________________________________
185 
186 void StFtpcLaserCalib::minuit_print()
187 {
188  // Print results
189  Double_t amin,edm,errdef;
190  Int_t nvpar,nparx,icstat;
191  LOG_INFO<<"StFtpcLaserCalib::minuit_print():"<<endm;
192  mMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
193  mMinuit->mnprin(5,amin);
194  LOG_INFO<<" "<<endm;
195 }
196 
197 //______________________________________________________________________________
198 
199 
200 void StFtpcLaserCalib::calc_res()
201 {
202  //LOG_DEBUG<<"Calculate residuals in x,y & r,phi ..."<<endm;
203  //LOG_DEBUG<<" "<<endm;
204 
205 
206  TF1 *fxz=new TF1("fxz","[0]*x+[1]",0,6);
207  TF1 *fyz=new TF1("fyz","[0]*x+[1]",0,6);
208 
209  // set fitted parameter to lin. function
210  Double_t v,ev;
211  for (int i=0;i<4;i++)
212  {
213  mMinuit->GetParameter(i,v,ev);
214  if (i<2) fxz->SetParameter(i,v);
215  else fyz->SetParameter(i-2,v);
216  }
217 
218  // calc residuals in x and y
219 
220  hnhits->Fill(nhits);
221 #ifdef HELIX_FIT
222  htrcharge->Fill(trcharge);
223  hp->Fill(p);
224  hpt->Fill(pt);
225  hinvp->Fill(1/p*trcharge);
226  hinvpt->Fill(1/pt*trcharge);
227 #endif
228 
229  for (int i=0;i<nhits;i++)
230  {
231  calcx[i]=fxz->Eval(z[i]);
232  calcy[i]=fyz->Eval(z[i]);
233 
234  usedfit++;
235 
236  resx[i]=(x[i]-fxz->Eval(z[i]));
237  resy[i]=(y[i]-fyz->Eval(z[i]));
238  calcrad[i]=sqrt(fxz->Eval(z[i])*fxz->Eval(z[i])+fyz->Eval(z[i])*fyz->Eval(z[i]));
239  if (FTPC==1)
240  rpol[i]=sqrt(fxz->Eval(zfieldcage)*fxz->Eval(zfieldcage)+fyz->Eval(zfieldcage)*fyz->Eval(zfieldcage));
241  if (FTPC==2)
242  rpol[i]=sqrt(fxz->Eval(-zfieldcage)*fxz->Eval(-zfieldcage)+fyz->Eval(-zfieldcage)*fyz->Eval(-zfieldcage));
243  calcphi[i]=zyltrafo(fxz->Eval(z[i]),fyz->Eval(z[i]),z[i]);
244  resrad[i]=radius[i]-calcrad[i];
245  resphi[i]=(phi[i]-calcphi[i])*Rad2Grad;
246 
247  }
248 
249  delete fxz;
250  delete fyz;
251 }
252 
253 //______________________________________________________________________________
254 
255 
256 int StFtpcLaserCalib::laser_fit(int getnhits)
257 {
258 
259  nhits=getnhits;
260 
261  int new_nhits=0;
262 
263  for (int i=0;i<nhits;i++)
264  {
265  if (cut(i))
266  {
267  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]);
268  new_nhits++;
269  }
270  }
271 
272  nhits=new_nhits;
273  if (nhits <= 2 ) {
274  //LOG_WARN << "StFtpcLaserCalib::laser_fit - Can't fit a laser track with " << nhits << " hits" << endm;
275  return -1;
276  }
277 
278  // ***************************************************
279  // * calculate new space point depending on gas & t0 *
280  // ***************************************************
281 
282 
283  //if ((deltat0!=0 || deltagas!=0) && nhits>0)
284  if (nhits>0)
285  {
286  //LOG_DEBUG<<"calculating new space point..."<<endm;
287  for (int i=0;i<nhits;i++) // debug : <=nhits
288  {
289 
290  //LOG_DEBUG<<sqrt(x[i]*x[i]+y[i]*y[i])<<endm;
291  //LOG_DEBUG<<softrow[i]<<" "<<softsec[i]<<" "<<timepos[i]<<" "<<ppos[i]<<" "<<x[i]<<" "<<y[i]<<endm;
292 
293  if (mtrafo->padtrans(softrow[i],softsec[i],timepos[i],ppos[i],&x[i],&y[i]))
294  {
295 
296  radius[i]=sqrt(x[i]*x[i]+y[i]*y[i]);
297  //LOG_DEBUG<<"radius["<<i<<"] = "<<radius[i]<<endm;
298  phi[i]=zyltrafo(x[i],y[i],z[i]);
299  }
300  else {
301  LOG_ERROR<<"ERROR padtrans !"<<endm;
302  }
303  }
304  //LOG_DEBUG<<"after new space point calculation !"<<endm;
305  }
306 
307 #ifdef HELIX_FIT
308  // ************************************************************************
309  // * include Helix fit & B-field distortion corrections (StFtpcTrack) !!!) *
310  // * keine Trafo noetig da nur im Global, ggf. fuer East !??? *
311  // ************************************************************************
312 
313  // create Laser Dummy Vertex (fuer data aus Tree auslesen !!! oder global nur);
314 
315  StFtpcVertex *laser_vertex=new StFtpcVertex(0,0,0,0,0,0);
316 
317  // create new track
318 
319  StFtpcTrack *mLaserTrack=new StFtpcTrack(1);
320 
321  TObjArray *mTrackHits=new TObjArray(0);
322 
323  // create new point and add to track
324  for (int i=0;i<nhits;i++)
325  {
326 // 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);
327  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);
328  mTrackHits->AddLast(mNewPoint);
329  //LOG_DEBUG<<mNewPoint->GetX()<<" "<<mNewPoint->GetY()<<" "<<mNewPoint->GetZ()<<endm;
330  mLaserTrack->AddPoint((StFtpcPoint*) mTrackHits->Last());
331  //delete mNewPoint;
332  }
333 
334  // DEBUG :
335  //LOG_DEBUG<<mLaserTrack->GetTrackNumber()<<" "<<mLaserTrack->GetNumberOfPoints()<<endm;
336  //LOG_DEBUG<<laser_vertex->GetX()<<" "<<laser_vertex->GetY()<<" "<<laser_vertex->GetZ()<<endm;
337 
338  // ********************************
339  // * Fit laser track with Helix ! *
340  // ********************************
341 
342 // mLaserTrack->Fit(laser_vertex,10.,true);
343  mLaserTrack->Fit(laser_vertex,10.,false);
344  if (mLaserTrack->GetP() == 0 && mLaserTrack->GetPt() == 0) {
345  //LOG_WARN << "StFtpcLaserCalib::laser_fit - Fit laser track with Helix was unsuccessful" << endm;
346  return -1;
347  }
348 
349  // DEBUG :
350  //LOG_DEBUG<<"p = "<<mLaserTrack->GetP()<<" | pt = "<<mLaserTrack->GetPt()<<" | charge = "<<mLaserTrack->GetCharge()<<endm;
351 
352  // get momentum and charge sign of Helix Fit !!!
353  trcharge=mLaserTrack->GetCharge();
354  p=mLaserTrack->GetP();
355  pt=mLaserTrack->GetPt();
356  //pt=1/(mLaserTrack->curvature());
357  //LOG_DEBUG<<mLaserTrack->curvature()<<endm;
358 
359  // delete Hit Array and get residuals before ...
360  for (int i=0;i<nhits;i++)
361  {
362 
363  resx2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetXGlobResidual();
364  resy2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetYGlobResidual();
365  resrad2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetRGlobResidual();
366  resphi2[i]=((StFtpcPoint*) mTrackHits->At(i))->GetPhiGlobResidual();
367 
368  // DEBUG :
369  //LOG_DEBUG<<" i = "<<i<<endm;
370  //LOG_DEBUG<<((StFtpcPoint*) mTrackHits->At(i))->GetXS()<<" "<<((StFtpcPoint*) mTrackHits->At(i))->GetX()<<endm;
371  //LOG_DEBUG<<((StFtpcPoint*) mTrackHits->At(i))->GetZS()<<" "<<((StFtpcPoint*) mTrackHits->At(i))->GetZ()<<endm;
372  // noch nicht alle Punkte richtig !!!!
373  // letzter Punkt wird nicht gesetzt (kleinstes z) vgl StFtpcTrack.cc !!!
374  // nehme dann standard !!!
375 
376  //LOG_DEBUG<<x[i]-((StFtpcPoint*) mTrackHits->At(i))->GetXS()<<endm;
377 
378  // create shiftet space position !!!
379 
380 // if (i<nhits-1)
381 // {
382 // x_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetXS();
383 // y_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetYS();
384 // }
385 // else
386 // {
387  x_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetX();
388  y_s[i]=((StFtpcPoint*) mTrackHits->At(i))->GetY();
389 // }
390 
391  //LOG_DEBUG<<x_s[i]-x[i]<<endm;
392 
393  // calc difference as check !!!
394  x_d[i]=x[i]-x_s[i];y_d[i]=y[i]-y_s[i];
395 
396  delete (StFtpcPoint*) mTrackHits->At(i);
397  }
398 
399  delete mTrackHits;
400 
401  delete laser_vertex;
402  delete mLaserTrack;
403 
404  // ***********************************
405  // * Straight line Fit with Miniut ! *
406  // ***********************************
407 
408 
409  // if gas or t0 changes take B-field corrected space pos.
410  // (check if chain default = corrected !?)
411 
412  //if ((deltat0!=0 || deltagas!=0) && nhits>0)
413 
414 
415  for (int i=0;i<nhits;i++)
416  {
417  x[i]=x_s[i];
418  y[i]=y_s[i];
419  }
420 
421 #endif
422 
423  minuit_init();
424  minuit_set_par();
425  //mMinuit->SetPrintLevel(-1);
426 
427  arglist[0] = 500;
428  arglist[1] = 1;
429 
430  if (nhits>0)
431  {
432  //LOG_DEBUG<<"start fit"<<endm;
433  mMinuit->mnexcm("MIGRAD", arglist , 2 ,ierflg);
434  //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm;
435  //check if fit okay !
436  if (ierflg==0)
437  {
438  //minuit_print();
439  calc_res();
440 
441  // Fill histograms !!!!
442  for (int i=0;i<nhits;i++)
443  Fill(i);
444 
445  }
446  }
447  else
448  {
449  LOG_ERROR<<"StFtpcLaserCalib::laser_fit - Number cluster on tracks = 0 !"<<endm;
450  ierflg=-1;
451  }
452 
453  //LOG_DEBUG<<"Minuit Error-Code : "<<ierflg<<endm;
454 
455  nhits=0;
456  return ierflg;
457 }
458 
459 //______________________________________________________________________________
460 
461 // erstmal test fuer Spectrum !
462 void StFtpcLaserCalib::defl_angle_transv()
463 {
464  TSpectrum *specrad=new TSpectrum();
465  //LOG_DEBUG<<"*******************************"<<endm;
466  //LOG_DEBUG<<"Anzahl der Hits (rad) : "<<specrad->Search(hrad,3)<<endm;
467  LOG_INFO<<"*******************************"<<endm;
468  for (int i=0;i<specrad->Search(hrad,2.5);i++)
469  {
470  LOG_INFO<<"Radius "<<i<<" : "<<specrad->GetPositionX()[i]<<endm;
471  }
472  LOG_INFO<<"*******************************"<<endm;
473  specrad->Delete();
474 }
475 
476 //______________________________________________________________________________
477 
478 void StFtpcLaserCalib::defl_histograms_st()
479 {
480  char name[30];
481  char name2[30];
482  for (int i=1;i<11;i++)
483  {
484  sprintf(name,"rad_l%d",i);
485  sprintf(name2,"pad_l%d",i);
486  if (i<4)
487  {
488  hpadcut[i]=new TH1F(name2,name2,200,105,125);
489  hradcut[i]=new TH1F(name,name,(int)((radcutmaxst[i]-radcutminst[i])/0.1),radcutminst[i],radcutmaxst[i]);
490  }
491  else
492  {
493  hpadcut[i]=new TH1F();
494  hradcut[i]=new TH1F();
495  }
496  }
497 }
498 
499 //______________________________________________________________________________
500 
501 void StFtpcLaserCalib::defl_histograms()
502 {
503  char name[30];
504  char name2[30];
505  for (int i=1;i<11;i++)
506  {
507  sprintf(name,"rad_l%d",i);
508  sprintf(name2,"pad_l%d",i);
509 
510  hpadcut[i]=new TH1F(name2,name2,300,120,150);
511  hradcut[i]=new TH1F(name,name,(int)((radcutmax[i]-radcutmin[i])/0.1),radcutmin[i],radcutmax[i]);
512  }
513 }
514 
515 //______________________________________________________________________________
516 
517 
518 void StFtpcLaserCalib::extrapol_histograms()
519 {
520  char name[30];
521 
522  for (int i=1;i<4;i++)
523  {
524  sprintf(name,"radpol_l%d",i);
525  hradpolcut[i]=new TH1F(name,name,(int)((radcutmaxst[i]-radcutminst[i])/0.025),radcutminst[i],radcutmaxst[i]);
526  }
527 }
528 
529 //______________________________________________________________________________
530 
531 void StFtpcLaserCalib::fill_extrapol_histograms(float getradpol)
532 {
533  for (int i=1;i<4;i++)
534  {
535  if (getradpol>radcutminst[i] && getradpol<radcutmaxst[i])
536  {
537  ((TH1F*) hradpolcut[i])->Fill(getradpol);
538  }
539  }
540 }
541 
542 //______________________________________________________________________________
543 
544 void StFtpcLaserCalib::fill_defl_histograms(float getrad, float getpadpos)
545 {
546  for (int i=1;i<11;i++)
547  {
548  if (getrad>radcutmin[i] && getrad<radcutmax[i])
549  {
550  ((TH1F*) hradcut[i])->Fill(getrad);
551  ((TH1F*) hpadcut[i])->Fill(getpadpos);
552 
553  }
554  }
555 }
556 
557 //______________________________________________________________________________
558 
559 
560 void StFtpcLaserCalib::fill_defl_histograms_st(float getrad, float getpadpos)
561 {
562  for (int i=1;i<4;i++)
563  {
564  if (getrad>radcutminst[i] && getrad<radcutmaxst[i])
565  {
566  ((TH1F*) hradcut[i])->Fill(getrad);
567  ((TH1F*) hpadcut[i])->Fill(getpadpos);
568 
569  }
570  }
571 }
572 
573 //______________________________________________________________________________
574 
575 void StFtpcLaserCalib::analyse_defl()
576 {
577  ofstream padfile,padfile2;
578  padfile.open(filename+"_padpos.log",ios::out);
579  padfile2.open(filename+"_padpos2.log",ios::out);
580 
581  // use mean
582  padfile<<"Padposition and radius with mean of histograms :"<<endl;
583  padfile<<"------------------------------------------------"<<endl;
584  padfile<<endl;
585  padfile<<"radius [cm] | padposition"<<endl;
586  padfile<<"-------------------------"<<endl;
587  padfile<<endl;
588 
589  for (int i=1;i<11;i++)
590  {
591  if (((TH1F*) hradcut[i])->GetEntries()>0)
592  {
593  padfile<<((TH1F*) hradcut[i])->GetMean()<<"+-"<<((TH1F*) hradcut[i])->GetRMS()<<" ";
594  padfile<<((TH1F*) hpadcut[i])->GetMean()<<"+-"<<((TH1F*) hpadcut[i])->GetRMS()<<endl;
595  padfile2<<((TH1F*) hradcut[i])->GetMean()<<" "<<((TH1F*) hradcut[i])->GetRMS()<<" ";
596  padfile2<<((TH1F*) hpadcut[i])->GetMean()<<" "<<((TH1F*) hpadcut[i])->GetRMS()<<endl;
597  }
598  }
599 
600  padfile<<endl;
601 
602  padfile<<"Padposition and radius with mean of gausfit :"<<endl;
603  padfile<<"------------------------------------------------"<<endl;
604  padfile<<endl;
605  padfile<<"radius [cm] | padposition"<<endl;
606  padfile<<"-------------------------"<<endl;
607  padfile<<endl;
608 
609  // use gausfit !!!
610  TF1 *gausfit3=new TF1("gausfit3","gaus");
611 
612  for (int i=1;i<11;i++)
613  {
614  if (((TH1F*) hradcut[i])->GetEntries()>0)
615  {
616  ((TH1F*) hradcut[i])->Fit("gausfit3","NQ");
617  padfile<<gausfit3->GetParameter(1)<<"+-"<<gausfit3->GetParameter(2)<<" ";
618  padfile2<<gausfit3->GetParameter(1)<<" "<<gausfit3->GetParameter(2)<<" ";
619  ((TH1F*) hpadcut[i])->Fit("gausfit3","NQ");
620  padfile<<gausfit3->GetParameter(1)<<"+-"<<gausfit3->GetParameter(2)<<endl;
621  padfile2<<gausfit3->GetParameter(1)<<" "<<gausfit3->GetParameter(2)<<endl;
622  }
623  }
624 
625  // make ps file !!!
626 
627 
628  padfile.close();
629  padfile2.close();
630 }
631 
632 //______________________________________________________________________________
633 
634 void StFtpcLaserCalib::PositionLog()
635 {
636 
637  ofstream logfile;
638  logfile.open(filename+".log",ios::out);
639 
640  TSpectrum *spec=new TSpectrum();
641  //LOG_DEBUG<<" "<<endm;
642  //LOG_DEBUG<<"Gaus fit to radius distribution of the "<<spec->Search(hrad,6)<<" reconstructed straight Laser tracks..."<<endm;
643  //LOG_DEBUG<<"-----------------------------------------------------------------------------"<<endm;
644  for (int i=0;i<spec->Search(hrad,6);i++)
645  {
646  //LOG_DEBUG<<"Radius after Peakfinder : "<<(spec->GetPositionX())[i]<<endm;
647  logfile<<"Radius after Peakfinder : "<<(spec->GetPositionX())[i]<<endl;
648  if (GAUSFIT==1)
649  {
650  float min=(spec->GetPositionX()[i]-1);
651  float max=(spec->GetPositionX()[i]+1);
652  TF1 *gausfit=new TF1("gausfit","gaus",min,max);
653  gausfit->SetParameter(1,(spec->GetPositionX())[i]);
654  // ggf. ParLimits einfuehren um richtige Fehlerbehandlung !!!
655  hrad->Fit(gausfit,"RQN");
656  //LOG_DEBUG<<"Radius after Gausfit : "<<gausfit->GetParameter(1)<<" +- "<<gausfit->GetParameter(2)<<endm;
657  logfile<<"Radius after Gausfit : "<<gausfit->GetParameter(1)<<" +- "<<gausfit->GetParameter(2)<<endl;
658  }
659  }
660  //LOG_DEBUG<<" "<<endm;
661  //LOG_DEBUG<<"Gaus fit to phi distribution of the "<<spec->Search(hphi,2)<<" reconstructed straight Laser tracks..."<<endm;
662  //LOG_DEBUG<<"-----------------------------------------------------------------------------"<<endm;
663  for (int i=0;i<spec->Search(hphi,6);i++)
664  {
665  //LOG_DEBUG<<"Phi after Peakfinder : "<<(spec->GetPositionX())[i]<<endm;
666  logfile<<"Phi after Peakfinder : "<<(spec->GetPositionX())[i]<<endl;
667  if (GAUSFIT==1)
668  {
669  float min=(spec->GetPositionX()[i]-0.5);
670  float max=(spec->GetPositionX()[i]+0.5);
671  TF1 *gausfit2=new TF1("gausfit2","gaus",min,max);
672  gausfit2->SetParameter(1,(spec->GetPositionX())[i]);
673  hphi->Fit(gausfit2,"RQN");
674  //LOG_DEBUG<<"Phi after Gausfit : "<<gausfit2->GetParameter(1)<<" +- "<<gausfit2->GetParameter(2)<<endm;
675  logfile<<"Phi after Gausfit : "<<gausfit2->GetParameter(1)<<" +- "<<gausfit2->GetParameter(2)<<endl;
676  }
677  }
678  logfile.close();
679 }
680 
681 //______________________________________________________________________________
682 
683 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)
684 {
685 
686  radius[n]=sqrt(tx*tx+ty*ty);
687  phi[n]=zyltrafo(tx,ty,tz);
688  //phi[n]=atan2(ty,tx);
689  x[n]=tx;
690  y[n]=ty;
691  z[n]=tz;
692  hsec[n]=nsec;
693  ex[n]=tex;
694  ey[n]=tey;
695  ppos[n]=gppos;
696  ppossigma[n]=gppossigma;
697  softsec[n]=gsoftsec;
698  softrow[n]=gsoftrow;
699  timepos[n]=gtimepos;
700  timel[n]=gettimel;
701  padl[n]=getpadl;
702  charge[n]=getcharge;
703  maxadc[n]=getmaxadc;
704 
705  //LOG_DEBUG<<radius[n]<<" "<<maxadc[n]<<endm;
706 }
707 
708 //______________________________________________________________________________
709 
710 void StFtpcLaserCalib::Fill(int l)
711 {
712  hrad->Fill(radius[l]);
713  hradpol->Fill(rpol[l]);
714  hphi->Fill(phi[l]*rad2grad);
715  hcalcrad->Fill(calcrad[l]);
716  hcalcphi->Fill(calcphi[l]);
717  hphiz->Fill(z[l],phi[l]);
718  hradz->Fill(z[l],radius[l]);
719  hresx->Fill(resx[l]);
720  hresy->Fill(resy[l]);
721  hresrad->Fill(resrad[l]);
722  hresphi->Fill(resphi[l]);
723  hresrad2->Fill(resrad[l],radius[l]);
724  hresphi2->Fill(resphi[l],radius[l]);
725  hresx2->Fill(resx[l],radius[l]);
726  hresy2->Fill(resy[l],radius[l]);
727 
728 #ifdef HELIX_FIT
729  hhresx->Fill(resx2[l]);
730  hhresy->Fill(resy2[l]);
731  hhresrad->Fill(resrad2[l]);
732  hhresphi->Fill(resphi2[l]*rad2grad);
733  hhresrad2->Fill(resrad2[l],radius[l]); // radius ggf. da Helix mit schieben anders !????
734  hhresphi2->Fill(resphi2[l]*rad2grad,radius[l]);
735  hhresx2->Fill(resx2[l],radius[l]);
736  hhresy2->Fill(resy2[l],radius[l]);
737 #endif
738 
739  hpad->Fill(ppos[l]);
740  htime->Fill(timepos[l]);
741  hpadrad->Fill(radius[l],ppos[l]);
742  hpadsigma->Fill(ppossigma[l]);
743  hpadl->Fill(padl[l]);
744  htimel->Fill(timel[l]);
745  hcharge->Fill(charge[l]);
746  hmaxadc->Fill(maxadc[l]);
747 
748 #ifdef HELIX_FIT
749  hbdiffx->Fill(x_d[l]);
750  hbdiffy->Fill(y_d[l]);
751 #endif
752 
753  if (STRAIGHT==0 || STRAIGHT==3)
754  fill_defl_histograms(radius[l],ppos[l]);
755  else if (STRAIGHT==1)
756  fill_defl_histograms_st(radius[l],ppos[l]);
757 
758  fill_extrapol_histograms(rpol[l]);
759 }
760 
761 
762 //______________________________________________________________________________
763 
764 void StFtpcLaserCalib::MakePs()
765 {
766  LOG_INFO<<"Create "<<filename+".ps"<<endm;
767 
768  gStyle->SetPalette(1);
769 
770  TCanvas *c1 = new TCanvas("c1","ps",200,10,700,500);
771  TPostScript *fps=new TPostScript(filename+".ps",112);
772 
773 
774  fps->NewPage();
775  c1->Divide(2,2);
776  c1->Update();
777  c1->cd(1);
778  hrad->DrawCopy();
779  c1->cd(2);
780  //hphi->DrawCopy();
781  hnhits->DrawCopy();
782  c1->cd(3);
783  hradz->DrawCopy("box");
784  c1->cd(4);
785  hphi->DrawCopy();
786  //hphiz->DrawCopy("box");
787  c1->Update();
788  fps->NewPage();
789  c1->cd(1);
790 
791  TF1 *gfit=new TF1("gfit","gaus",-0.5,0.5);
792  gfit->FixParameter(1,0.0);
793  //gfit->SetParameter(1,0.0);
794  //gfit->SetParLimits(1,0.0,0.0);
795  gfit->SetLineColor(2);
796  gfit->SetLineWidth(1);
797 
798  ofstream logfile;
799  logfile.open(filename+"_res.log",ios::out);
800  c1->cd(1);
801  hresx->DrawCopy();
802  hresx->Fit(gfit,"QR");
803  logfile<<deltat0<<" "<<deltagas<<" "<<deltaTemp<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl;
804  c1->cd(2);
805  hresy->DrawCopy();
806  hresy->Fit(gfit,"QR");
807  logfile<<deltat0<<" "<<deltagas<<" "<<deltaTemp<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl;
808  c1->cd(3);
809  hresrad->DrawCopy();
810  hresrad->Fit(gfit,"QR");
811  logfile<<deltat0<<" "<<deltagas<<" "<<deltaTemp<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl;
812  c1->cd(4);
813  hresphi->DrawCopy();
814  hresphi->Fit(gfit,"QR");
815  logfile<<deltat0<<" "<<deltagas<<" "<<deltaTemp<<" "<<gfit->GetParameter(2)<<" "<<gfit->GetChisquare()<<endl;
816  //logfile<<"Delta t0 | Delta Gas | Sigma res (x,y,r,phi)"<<endl;
817  //logfile.close();
818 
819  c1->Update();
820  fps->NewPage();
821  c1->cd(1);
822  hresrad2->DrawCopy("colz");
823  c1->cd(2);
824  hresphi2->DrawCopy("colz");
825  //c1->cd(3);
826  //hrad->DrawCopy();hcalcrad->DrawCopy("same");
827  //c1->cd(4);
828  //hphi->DrawCopy();hcalcphi->DrawCopy("same");
829  c1->cd(3);
830  hresx2->DrawCopy("colz");
831  c1->cd(4);
832  hresy2->DrawCopy("colz");
833  c1->Update();
834  fps->NewPage();
835  c1->cd(1);
836  hradpol->Draw();
837  //logfile<<deltat0<<" "<<deltagas<<" "<<hradpol->GetMean()<<endl;
838  //logfile.close();
839  c1->cd(2);
840  hpad->DrawCopy();
841  c1->cd(3);
842  hpadrad->DrawCopy("colz");
843  c1->cd(4);
844  hpadsigma->DrawCopy();
845  c1->Update();
846  fps->NewPage();
847 
848  for (int i=1;i<4;i++)
849  {
850  c1->cd(i);
851  ((TH1F*) hradpolcut[i])->DrawCopy();
852  logfile<<deltat0<<" "<<deltagas<<" "<<deltaTemp<<" "<<((TH1F*) hradpolcut[i])->GetMean()<<" "<<((TH1F*) hradpolcut[i])->GetRMS()<<endl;
853  logfile<<deltat0<<" "<<deltagas<<" "<<deltaTemp<<" "<<(((TH1F*) hradpolcut[i])->GetXaxis())->GetBinCenter(((TH1F*) hradpolcut[i])->GetMaximumBin())<<endl;
854  }
855  c1->Update();
856 
857  fps->NewPage();
858  c1->cd(1);hpadl->DrawCopy();
859  c1->cd(2);htimel->DrawCopy();
860  c1->cd(3);hmaxadc->DrawCopy();//hmaxadc->Fit("landau");
861  c1->cd(4);hcharge->DrawCopy();//hcharge->Fit("landau");
862  c1->Update();
863 
864 #ifdef HELIX_FIT
865  fps->NewPage();
866  //c1->cd(1);htrcharge->DrawCopy();
867  c1->cd(1);hinvp->DrawCopy();hinvp->Fit(gfit,"QR");
868  c1->cd(2);hinvpt->DrawCopy();hinvpt->Fit(gfit,"QR");
869  c1->cd(3);hbdiffx->DrawCopy();
870  c1->cd(4);hbdiffy->DrawCopy();
871  c1->Update();
872 
873  gfit->SetLineColor(3);
874 
875  fps->NewPage();
876  c1->cd(1);
877  hhresx->DrawCopy();
878  hhresx->Fit(gfit,"QR");
879  c1->cd(2);
880  hhresy->DrawCopy();
881  hhresy->Fit(gfit,"QR");
882  c1->cd(3);
883  hhresrad->DrawCopy();
884  hhresrad->Fit(gfit,"QR");
885  c1->cd(4);
886  hhresphi->DrawCopy();
887  hhresphi->Fit(gfit,"QR");
888  c1->Update();
889 
890  fps->NewPage();
891  c1->cd(1);
892  hhresrad2->DrawCopy("colz");
893  c1->cd(2);
894  hhresphi2->DrawCopy("colz");
895  c1->cd(3);
896  hhresx2->DrawCopy("colz");
897  c1->cd(4);
898  hhresy2->DrawCopy("colz");
899  c1->Update();
900 
901  fps->NewPage();
902  c1->cd(1);
903  hhresx->DrawCopy();
904  hresx->SetLineColor(2);
905  hresx->DrawCopy("same");
906  c1->cd(2);
907  hhresy->DrawCopy();
908  hresy->SetLineColor(2);
909  hresy->DrawCopy("same");
910  c1->cd(3);
911  hhresrad->DrawCopy();
912  hresrad->SetLineColor(2);
913  hresrad->DrawCopy("same");
914  c1->cd(4);
915  hhresphi->DrawCopy();
916  hresphi->SetLineColor(2);
917  hresphi->DrawCopy("same");
918  c1->Update();
919 #endif
920 
921  logfile.close();
922  fps->Close();
923  //c1->Delete();
924 }
925 
926 //______________________________________________________________________________
927 
928 void StFtpcLaserCalib::MakeOutput(TString eingabe,char* t0, char* gas, float gastemp)
929 {
930 
931 stringstream dT;
932 dT << gastemp;
933 std::string dTstring = dT.str();
934 const char *dTemp = dTstring.c_str();
935 
936  eingabe +="-pos_";
937  if (FTPC==1)
938  eingabe +="w_lsec_";
939  else if (FTPC==2)
940  eingabe +="e_lsec_";
941  else if (FTPC==0)
942  eingabe +="all_lsec";
943 
944  if (FTPC!=0)
945  eingabe +=LSEC;
946 
947  if (STRAIGHT==1)
948  eingabe+="_g";
949  else if (STRAIGHT==0)
950  eingabe+="_s";
951  else if (STRAIGHT==2)
952  eingabe+="_s2";
953  else if (STRAIGHT==3)
954  eingabe+="_all";
955 
956  eingabe +="_dt"; eingabe +=t0;
957  eingabe +="_dg"; eingabe +=gas;
958  eingabe +="_dT"; eingabe +=dTemp;
959 
960  filename=eingabe;
961 
962  LOG_INFO<<"Create file : "<<filename+".root"<<endm;
963  outfile=new TFile(filename+".root","RECREATE");
964 
965  hrad=new TH1F("rad","radius (straight) laser tracks",124*8,0.5,31.5);
966  hradpol=new TH1F("radpol","radius laser tracks extrapoliert fieldcage",124*8,0.5,31.5);
967  //hphi=new TH1F("phi","phi (straight) laser tracks",48*8,-6,6);
968  hphi=new TH1F("phi","phi laser tracks",360*4,-90,270);
969  hcalcrad=new TH1F("clacrad","calc radius",124*8,0.5,31.5);
970  hcalcphi=new TH1F("calcphi","calc phi",48*8,-6,6);
971  hcalcrad->SetLineColor(2);hcalcphi->SetLineColor(2);
972 
973  if (FTPC==1)
974  {
975  hradz=new TH2F("radz","radius (straight) laser tracks",100,140,275,124*4,0.5,31.5);
976  hphiz=new TH2F("phiz","phi (straight) laser tracks",100,140,275,48*2,0,4);
977  }
978  else if (FTPC==2)
979  {
980  hradz=new TH2F("radz","radius (straight) laser tracks",100,-275,-140,124*4,0.5,31.5);
981  hphiz=new TH2F("phiz","phi (straight) laser tracks",100,-275,-140,48*2,0,4);
982  }
983  else if (FTPC==0)
984  {
985  hradz=new TH2F("radz","radius (straight) laser tracks",200,-275,275,124*4,0.5,31.5);
986  hphiz=new TH2F("phiz","phi (straight) laser tracks",200,-275,275,48*2,0,4);
987  }
988 
989  Int_t r_vs_z_bin=124/2;
990 
991  hresx=new TH1F("resx","Residual x",100,-0.5,0.5);
992  hresy=new TH1F("resy","Residual y",100,-0.5,0.5);
993  hresrad=new TH1F("resrad","Residual radius",100,-0.5,0.5);
994  hresphi=new TH1F("resphi","Residual phi",100,-0.5,0.5);
995  hresx2=new TH2F("resx2","Residual x vs.radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
996  hresx2->SetMinimum(0);
997  hresy2=new TH2F("resy2","Residual y vs. radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
998  hresy2->SetMinimum(0);
999  hresrad2=new TH2F("resrad2","Residual radius vs.radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
1000  hresrad2->SetMinimum(0);
1001  hresphi2=new TH2F("resphi2","Residual phi vs. radius",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
1002  hresphi2->SetMinimum(0);
1003  hpad=new TH1F("pados","Padposition",1600,0,160);
1004  htime=new TH1F("timepos","Timeposition",1800,0,180);
1005  hpadrad=new TH2F("radpados","Padposition vs. radius",r_vs_z_bin,0.5,31.5,1600,0,160);
1006  hpadsigma=new TH1F("padsigma","Sigma padposition",50,0,5);
1007 
1008  hmaxadc=new TH1F("maxadc","MaxAdc",150,0.5,150.5);
1009  hcharge=new TH1F("charge","Charge",300,0.5,1500.5);
1010  hpadl=new TH1F("padl","Padlength",15,0.5,15.5);
1011  htimel=new TH1F("timel","Timelength",20,0.5,20.5);
1012 
1013  hnhits=new TH1F("nhits","Number hits on laser track",8,3.5,11.5);
1014 
1015 #ifdef HELIX_FIT
1016  htrcharge=new TH1F("trcharge","Charge of Helix Fit",3,-1.5,1.5);
1017  hp=new TH1F("p","p Helix Fit",100, -0.5, 0.5);
1018  hpt=new TH1F("pt","pt Helix Fit",100, -0.5, 0.5);
1019  hinvp=new TH1F("invp","1/p * charge of Helix Fit",100, -0.5, 0.5);
1020  hinvpt=new TH1F("invpt","1/pt * charge of Helix Fit",50, -0.5, 0.5);
1021 
1022  hhresx=new TH1F("hresx","Residual x Helix Fit",100,-0.5,0.5);
1023  hhresy=new TH1F("hrresy","Residual y Helix Fit",100,-0.5,0.5);
1024  hhresrad=new TH1F("hrresrad","Residual radius Helix Fit",100,-0.5,0.5);
1025  hhresphi=new TH1F("hresphi","Residual phi Helix Fit",100,-0.5,0.5);
1026  hhresrad2=new TH2F("hresrad2","Residual radius vs.radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
1027  hhresphi2=new TH2F("hresphi2","Residual phi vs. radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
1028  hhresx2=new TH2F("hresx2","Residual x vs.radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
1029  hhresy2=new TH2F("hresy2","Residual y vs. radius Helix Fit",100,-0.5,0.5,r_vs_z_bin,0.5,31.5);
1030 
1031  hbdiffx=new TH1F("bdiffx","Differenz x and B-field corr. x",100,-0.01,0.01);
1032  hbdiffy=new TH1F("bdiffy","Differenz y and B-field corr. y",100,-0.01,0.01);
1033  // if you activate the following two histograms be sure to activate the corresponding delete statement
1034  //hbdiffx2=new TH2F("bdiffx2","Differenz x and B-field corr. x vs. z",);
1035  //hbdiffy2=new TH2F("bdiffy2","Differenz y and B-field corr. y vs. z",);
1036 #endif
1037 
1038  if (STRAIGHT==0)
1039  defl_histograms();
1040  else
1041  defl_histograms_st();
1042 
1043  //LOG_DEBUG<<"For extrapolhistos"<<endm;
1044  extrapol_histograms();
1045 
1046 }
1047 
1048 
1049 //______________________________________________________________________________
1050 
1051 StFtpcLaserCalib::~StFtpcLaserCalib()
1052 {
1053  //if (deltat0=0)
1054  //delete trafo;
1055 
1056  outfile->Write();
1057 
1058  delete mMinuit;
1059 
1060  //LOG_DEBUG<<"for normal histos"<<endm;
1061 
1062  delete hrad;delete hradpol;delete hphi; delete hcalcrad;
1063  delete hcalcphi; delete hradz; delete hphiz;
1064  delete hresx;delete hresy;delete hresphi; delete hresrad; delete hresrad2;
1065  delete hresphi2; delete hpad; delete hpadrad; delete hpadsigma; delete htime;
1066  delete hnhits;
1067 #ifdef HELIX_FIT
1068  delete hhresx;delete hhresy;delete hhresphi; delete hhresrad; delete hhresrad2; delete hhresphi2;
1069  delete htrcharge;
1070  delete hp; delete hpt; delete hinvp; delete hinvpt;
1071  delete hhresx2; delete hhresy2;
1072  delete hresx2; delete hresy2;
1073  delete hbdiffx; delete hbdiffy;
1074  //delete hbdiffx2;delete hbdiffy2;
1075 #endif
1076 
1077  if (STRAIGHT==0)
1078  {
1079  for (int i=1;i<4;i++)
1080  {delete hpadcut[i];delete hradcut[i];}
1081  }
1082  else
1083  {
1084  for (int i=1;i<11;i++)
1085  {delete hpadcut[i];delete hradcut[i];}
1086  }
1087  outfile->Close();
1088  //LOG_DEBUG<<"StFtpcLaserCalib() deconstructed !"<<endm;
1089 }
Generated by   doxygen 1.8.5