TPC/fastSimul/testEdgeFit.C

   1 /// \file testEdgeFit.C
   2 ///
   3 /// ~~~
   4 /// .L $ALICE_ROOT/TPC/fastSimul/testEdgeFit.C+
   5 /// testEdge(50000);
   6 /// ~~~
   7
   8 #include "TVectorD.h"
   9 #include "TMath.h"
  10 #include "TLinearFitter.h"
  11 #include "TRandom.h"
  12 #include "TH2F.h"
  13 #include "TFile.h"
  14 #include "TTree.h"
  15 #include "TLegend.h"
  16 #include "TPad.h"
  17
  18 #include "TTreeStream.h"
  19 //#include "TVectorD.h"
  20
  21
  22 void testEdge(Int_t ntest){
  23
  24   TTreeSRedirector *cstream = new TTreeSRedirector("testEdge.root");
  25   Double_t yth=1;
  26   Double_t ymc0=0;
  27   Double_t kx=0;
  28   Double_t sy=0.25;
  29   Int_t npoints=100;
  30   TVectorD vecX;
  31   TVectorD vecY0;
  32   TVectorD vecYC;
  33   TVectorD vecYT;
  34
  35   TVectorD vecU;
  36   TVectorD vecYTa;
  37   TVectorD vecYTb;
  38   TVectorD vecYEa;
  39   TVectorD vecYEb;
  40   //
  41   TVectorD vecYTM5;
  42   TVectorD vecYTk5;
  43   //
  44   TVectorD vecE;
  45   //
  46   //
  47   for (Int_t it=0; it<ntest;it++){
  48     //
  49     ymc0 = gRandom->Rndm()*1;
  50     kx = gRandom->Rndm()*1.;
  51     npoints = TMath::Nint(gRandom->Rndm()*50);
  52     vecX.ResizeTo(npoints);
  53     vecY0.ResizeTo(npoints);
  54     vecYC.ResizeTo(npoints);
  55     vecYT.ResizeTo(npoints);
  56     //
  57     vecYTa.ResizeTo(npoints);
  58     vecYTb.ResizeTo(npoints);
  59     vecYEa.ResizeTo(npoints);
  60     vecYEb.ResizeTo(npoints);
  61     vecU.ResizeTo(npoints);
  62     //
  63     TLinearFitter fitter(2,"pol1");
  64     for (Int_t ipoint=npoints-1;ipoint>0;ipoint--){
  65       Double_t x[10];
  66       //
  67       x[0]=(Double_t)(ipoint)*0.2;
  68       Double_t y0   = ymc0+Double_t(x[0])*kx;
  69       Double_t yc   = y0+gRandom->Gaus()*sy;
  70       //
  71       if (yc<yth)   yc=yth;
  72       vecX[ipoint]   = x[0];;
  73       vecY0[ipoint]  = y0;
  74       vecYC[ipoint]  = yc;
  75       vecU[ipoint]   = 0;
  76       vecYT[ipoint]  = -10;
  77       vecYTa[ipoint] = -10;
  78       vecYTb[ipoint] = -10;
  79       fitter.AddPoint(x,yc,sy);
  80       if (fitter.GetNpoints()>4){
  81         fitter.Eval();
  82         fitter.GetErrors(vecE);
  83         vecYTa[ipoint]=fitter.GetParameter(0);
  84         vecYTb[ipoint]=fitter.GetParameter(1);
  85         vecYEa[ipoint]=vecE[0];
  86         vecYEb[ipoint]=vecE[1];
  87         vecU[ipoint]=fitter.GetNpoints();
  88         vecYT[ipoint]=fitter.GetParameter(0)+fitter.GetParameter(1)*x[0];
  89       }
  90     }
  91     (*cstream)<<"Dump"<<
  92       "n="<<npoints<<
  93       "ymc="<<ymc0<<
  94       "kx="<<kx<<
  95       "vX.="<<&vecX<<
  96       "vU.="<<&vecU<<
  97       "vY0.="<<&vecY0<<
  98       "vYC.="<<&vecYC<<
  99       "vYT.="<<&vecYT<<
 100       "vYTax.="<<&vecYTa<<
 101       "vYTb.="<<&vecYTb<<
 102       "vYEa.="<<&vecYEa<<
 103       "vYEb.="<<&vecYEb<<
 104       "\n";
 105   }
 106   delete cstream;
 107 }
 108
 109
 110 void MakePicDy(){
 111   TFile f("testEdge.root");
 112   TTree * tree = (TTree*)f.Get("Dump");
 113   TObjArray arrfdy0,arrfdyT;
 114   TH2F *phisdY0 = new TH2F("hisdY0","hisdY0",25,1,3,200,-3,3);
 115   TH2F *phisdYT = new TH2F("hisdYT","hisdYT",25,1,3,200,-3,3);
 116   tree->Draw("10*(vYT.fElements-vY0.fElements):vY0.fElements>>hisdY0","vU.fElements>20&&kx>0.1","colz");
 117   phisdY0->SetDirectory(0);
 118   tree->Draw("10*(vYT.fElements-vY0.fElements):vYT.fElements>>hisdYT","vU.fElements>20&&kx>0.1","colz");
 119   phisdYT->SetDirectory(0);
 120
 121   phisdY0->FitSlicesY(0,0,-1,0,"QNR",&arrfdy0);
 122   phisdYT->FitSlicesY(0,0,-1,0,"QNR",&arrfdyT);
 123   //
 124   TH1 *hisdy0 = (TH1*)arrfdy0.At(1);
 125   TH1 *hisdyT = (TH1*)arrfdyT.At(1);
 126   hisdy0->SetLineColor(1);
 127   hisdyT->SetLineColor(2);
 128   hisdy0->SetMarkerColor(1);
 129   hisdyT->SetMarkerColor(2);
 130   hisdy0->SetMarkerStyle(22);
 131   hisdyT->SetMarkerStyle(24);
 132   hisdyT->SetXTitle("Y (cm)");
 133   hisdyT->SetYTitle("Y_{t}-Y_{mc} (mm)");
 134   hisdyT->SetMinimum(-hisdyT->GetMaximum());
 135   hisdyT->Draw();
 136   hisdy0->Draw("same");
 137   TLegend *legend = new TLegend(0.25,0.60,0.85,0.85, "Track residuals close to edge (Edge= 1 cm, #sigma_{cl}=0.25 cm)");
 138   legend->AddEntry(hisdy0,"Y=Y_{mc}");
 139   legend->AddEntry(hisdyT,"Y=Y_{t}");
 140   legend->Draw();
 141   gPad->SaveAs("picEdge/dYedgeMC.eps");
 142   gPad->SaveAs("picEdge/dYedgeMC.gif");
 143 }
 144
 145 void MakePickDy(){
 146   TFile f("testEdge.root");
 147   TTree * tree = (TTree*)f.Get("Dump");
 148   TObjArray arrfdky0,arrfdkyT;
 149   TH2F *phisdKY0 = new TH2F("hisdKY0","hisdKY0",25,1,3,200,-150,150);
 150   TH2F *phisdKYT = new TH2F("hisdKYT","hisdKYT",25,1,3,200,-150,150);
 151   tree->Draw("1000*(vYTb.fElements-kx):vY0.fElements>>hisdKY0","vU.fElements>20","colz");
 152   phisdKY0->SetDirectory(0);
 153   tree->Draw("1000*(vYTb.fElements-kx):vYT.fElements>>hisdKYT","vU.fElements>20","colz");
 154   phisdKYT->SetDirectory(0);
 155
 156   phisdKY0->FitSlicesY(0,0,-1,0,"QNR",&arrfdky0);
 157   phisdKYT->FitSlicesY(0,0,-1,0,"QNR",&arrfdkyT);
 158   //
 159   TH1 *hisdky0 = (TH1*)arrfdky0.At(1);
 160   TH1 *hisdkyT = (TH1*)arrfdkyT.At(1);
 161   hisdky0->SetLineColor(1);
 162   hisdkyT->SetLineColor(2);
 163   hisdky0->SetMarkerColor(1);
 164   hisdkyT->SetMarkerColor(2);
 165   hisdky0->SetMarkerStyle(22);
 166   hisdkyT->SetMarkerStyle(24);
 167   hisdkyT->SetMinimum(-20);
 168   hisdkyT->SetMaximum(20);
 169   hisdkyT->SetXTitle("Y (cm)");
 170   hisdkyT->SetYTitle("k_{yt}-k_{ymc} (mrad)");
 171   hisdkyT->Draw();
 172   hisdky0->Draw("same");
 173   TLegend *legend = new TLegend(0.25,0.60,0.85,0.85, "Track residuals close to edge (Edge= 1 cm, #sigma_{cl}=0.25 cm)");
 174   legend->AddEntry(hisdky0,"Y=Y_{mc}");
 175   legend->AddEntry(hisdkyT,"Y=Y_{t}");
 176   legend->Draw();
 177   gPad->SaveAs("picEdge/dKYedgeMC.eps");
 178   gPad->SaveAs("picEdge/dKYedgeMC.gif");
 179 }
 180
 181