remove temporary task that was used to troubleshoot the gamma flow anaysis
[u/mrichter/AliRoot.git] / macros / TestRieman.C
CommitLineData
a45a4360 1#include "AliRieman.h"
2#include "TTreeStream.h"
3#include "TRandom.h"
4#include "AliExternalTrackParam.h"
5
6/*
7 //
8 // Test Progaram for AliRieman class
9 //
10 // How to use it:
11 // 1. Load compiled macros
12 // 2. Run function TestRieman()
13 // 3. Open file and tree with results
14 // 4. Check results - residuals - pulls
15 // See example bellow:
16
17 .L AliGenInfo.C+
18 .L TestRieman.C+
19 TestRieman();
20
21 TFile f("TestRieman.root");
22 TTree * tree = (TTree*)f.Get("Test");
23 tree->Draw("Rieman.fZ-RiemanRef.fZ:Rieman.fX")
24
25 AliComparisonDraw comp;
26 comp->fTree = tree;
27 comp->DrawXY("Par0.fP[2]","(Par0.fP[4]-ParR.fP[4])/sqrt(ParR.fC[14])","1","1",5,-0.1,0.1,-5.1,5.1);
28 comp->fRes->Draw();
29
30*/
31
32
33
34void GetProlongation(Double_t xk, Double_t &x, Double_t *param, Double_t &y, Double_t &z){
35 Double_t dx=xk-x;
36 Double_t f1=param[2], f2=f1 + param[4]*dx;
60e55aee 37 Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
a45a4360 38 y = param[0] + dx*(f1+f2)/(r1+r2);
39 //z = param[1] + dx*(r2 + f2*(f1+f2)/(r1+r2))*param[3];
40 //z = param[1] + dx*(f1+f2)/(f1*r2 + f2*r1)*param[3];
41 Double_t dy = y-param[0];
42 Double_t dfi = 2*TMath::ASin(0.5*TMath::Sqrt(dx*dx+dy*dy)*TMath::Abs(param[4]));
43 Double_t sign = (dx>0) ? 1:-1;
44 z = param[1] + sign*param[3]*dfi/param[4];
45}
46
47
48void TestRieman(){
49 const Double_t kB2C=0.299792458e-3;
50 TTreeSRedirector cstream("TestRieman.root");
51 AliRieman rieman(1000);
52 AliRieman riemanR(1000);
53 AliRieman riemanRef(1000);
54 Int_t npoints =150;
55 Double_t param[5], paramIdeal[5], paramR[5];
56
57 for (Int_t i=0; i<10000; i++){
58 //
59 // random samples
60 //
61 Double_t r = 600*(1+0.5*(gRandom->Rndm()-0.5));
62 Double_t kz = 2*(gRandom->Rndm()-0.5);
63 Double_t snp = (gRandom->Rndm()-0.5)*0.3;
64 Double_t sy =0.1, sz =0.1;
65 Double_t sign = (gRandom->Rndm()>0.5? -1:1);
66
67 Double_t x0 = 0;
68 Double_t y0 = snp*x0+10*(gRandom->Rndm()-0.5);
69 Double_t z0 = kz*x0+20*(gRandom->Rndm()-0.5);
70 param[0] = y0;
71 param[1] = z0;
72 param[2] = snp;
73 param[3] = kz;
74 param[4] = sign/r;
75 Double_t covar[15] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
76 Double_t covarI[15] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
77 Double_t covarR[15] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
78 AliExternalTrackParam trackParam(x0, 0, param, covar);
79 //
80 //
81 rieman.Reset();
82 riemanR.Reset();
83 riemanRef.Reset();
84 // for (Float_t dx =0; dx<90; dx+=5){trackParam.PropagateTo(dx,1./kB2C);}
85 for (Int_t ipoint =0; ipoint<npoints; ipoint++){
86 Double_t x = 90 + ipoint;
87 Double_t y,z;
88 GetProlongation(x, x0, param,y,z);
89 // trackParam.PropagateTo(x,1./kB2C);
90 // y = trackParam.GetY();
91 // z = trackParam.GetZ();
92 // for (Float_t dx =0; dx<1; dx+=1){trackParam.PropagateTo(x+dx,1./kB2C);}
93 // Double_t xyz[3];
94 // trackParam.GetXYZAt(x,1./kB2C,xyz);
95 // y = xyz[1];
96 // z = xyz[2];
97 //
98 rieman.AddPoint(x,y,z,sy,sz);
99 riemanR.AddPoint(x, gRandom->Gaus(y,sy), gRandom->Gaus(z,sz), sy,sz);
100 }
101 rieman.Update();
102 riemanR.Update();
103 //
104 // track extrapolation errrs
105 //
106 for (Int_t ipoint =0; ipoint<npoints; ipoint++){
107 Double_t x = rieman.GetX()[ipoint];
108 Double_t y = riemanR.GetYat(x);
109 Double_t z = riemanR.GetZat(x);
110 Double_t sty = riemanR.GetErrY(x);
111 Double_t stz = riemanR.GetErrZ(x);
112 riemanRef.AddPoint(x, y,z,sty,stz);
113 }
114 riemanRef.Update();
115 rieman.GetExternalParameters(x0, paramIdeal,covarI);
116 riemanR.GetExternalParameters(x0, paramR,covarR);
117
118 AliExternalTrackParam trackParamI(x0, 0, paramIdeal, covarI);
119 AliExternalTrackParam trackParamR(x0, 0, paramR, covarR);
120 AliRieman * res = rieman.MakeResiduals(); // residuals
121 AliRieman * resR = riemanR.MakeResiduals(); // residuals
122 cstream<<"Test"<<
123 "x0="<<x0<<
124 "Par0.="<<&trackParam<<
125 "ParI.="<<&trackParamI<<
126 "ParR.="<<&trackParamR<<
127 "Rieman.="<<&rieman<<
128 "RiemanR.="<<&riemanR<<
129 "RiemanRef.="<<&riemanRef<<
130 "Res.="<<res<<
131 "ResR.="<<resR<<
132 "\n";
133 delete res;
134 delete resR;
135 }
136}
137