macros/TestRieman.C

   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
  34 void 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;
  37   Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
  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
  48 void 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