TPC/fastSimul/testEdgeFit.C

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