ITS UPGRADE
[u/mrichter/AliRoot.git] / ITS / UPGRADE / ITSUpgradeRec / AliITSUMatLUT.cxx
1 #include "AliITSUMatLUT.h"
2 #include "AliLog.h"
3 #include "AliTrackerBase.h"
4 #include <TString.h>
5 #include <TH1.h>
6 #include <TMath.h>
7 #include <TRandom.h>
8
9
10 //___________________________________________________________
11 AliITSUMatLUT::AliITSUMatLUT() 
12   :fRMin(0)
13   ,fRMax(0)
14   ,fDRInv(-1)
15   ,fDR(0)
16   ,fNBins(0) 
17 {
18   // def c-tor
19   for (int i=kNParTypes;i--;) fData[i]=0;
20 }
21
22 //___________________________________________________________
23 AliITSUMatLUT::~AliITSUMatLUT() 
24 {
25   // d-tor
26   for (int i=kNParTypes;i--;) delete fData[i];
27 }
28
29 //___________________________________________________________
30 AliITSUMatLUT::AliITSUMatLUT(Double_t rmin,Double_t rmax,Int_t nbin)
31   :fRMin(rmin)
32   ,fRMax(rmax)
33   ,fDRInv(0)
34   ,fDR(0)
35   ,fNBins(nbin)
36 {
37   //
38   if (rmin<0 || rmax-rmin<1e-4 || nbin<1) AliFatal(Form("Illegal parameters Rmin:%f Rmax:%f Nbins:%d",rmin,rmax,nbin));
39   fDRInv = fNBins/(fRMax-fRMin);
40   fDR    = (fRMax-fRMin)/fNBins;
41   for (int i=kNParTypes;i--;) {
42     fData[i] = new Double_t[fNBins];
43     memset(fData[i],0,fNBins*sizeof(Double_t));
44   }
45   //
46 }
47
48 //___________________________________________________________
49 AliITSUMatLUT::AliITSUMatLUT(const AliITSUMatLUT& src)
50   :TObject(src)
51   ,fRMin(src.fRMin)
52   ,fRMax(src.fRMax)
53   ,fDRInv(src.fDRInv)
54   ,fDR(src.fDR)
55   ,fNBins(src.fNBins)
56 {
57   //
58   for (int i=kNParTypes;i--;) {
59     if (src.fData[i]) {
60       fData[i] = new Double_t[fNBins];
61       memcpy(fData[i],src.fData[i],fNBins*sizeof(Double_t));
62     }
63   }
64   //
65 }
66
67 //___________________________________________________________
68 AliITSUMatLUT & AliITSUMatLUT::operator=(const AliITSUMatLUT& src)
69 {
70   // copy 
71   if (this == &src) return *this;
72   this->~AliITSUMatLUT();
73   new(this) AliITSUMatLUT(src);
74   return *this;
75   //  
76 }
77
78 //___________________________________________________________
79 void AliITSUMatLUT::FillData(Int_t ntest, Double_t zmin,Double_t zmax)
80 {
81   // filla material data 
82   double start[3],stop[3],parStep[7];
83   if (fNBins<1) AliFatal("Limits are not set");
84   if (ntest<1 || zmin>zmax) AliFatal(Form("Wrong parameters Ntest:%d Zmin:%f Zmax:%f",ntest,zmin,zmax));
85   double dr = (fRMax-fRMin)/fNBins;
86   AliInfo(Form("Building material table for %.3f<R<%.3f %.3f<Z<%.3f in %d bins using %d tracks",fRMin,fRMax,zmin,zmax,fNBins,ntest));
87   const double kAngEps = 1e-4; // tiny slope to avoid tracks strictly normal to Z axis
88   double *tmpAcc = new double[fNBins*kNParTypes];
89   for (int itst=ntest;itst--;) {
90     double parInt[kNParTypes]={0};
91     double r   = fRMin;
92     double phi = gRandom->Rndm()*TMath::Pi()*2;
93     double cs  = TMath::Cos(phi);
94     double sn  = TMath::Sin(phi);
95     double angz = 2*(gRandom->Rndm()-0.5)*kAngEps;
96     stop[0] = r*cs;
97     stop[1] = r*sn;
98     stop[2] = zmin + gRandom->Rndm()*(zmax-zmin);
99     Bool_t fail = kFALSE;
100     for (int ir=0;ir<fNBins;ir++) {
101       r += dr;
102       for (int i=3;i--;) start[i] = stop[i];
103       stop[0] = r*cs;
104       stop[1] = r*sn;
105       stop[2] += dr*angz; // to avoid tracks normal to axis
106       AliTrackerBase::MeanMaterialBudget(start,stop, parStep);
107       if (parStep[1]>999) {fail = kTRUE; printf("fail\n"); break;}
108       //
109       parInt[kParX2X0] += parStep[1];
110       parInt[kParRhoL] += parStep[0]*parStep[4];
111       //
112       for (int ip=kNParTypes;ip--;) tmpAcc[ir*kNParTypes+ip] = parInt[ip];
113     }
114     if (fail) {itst++; continue;} // propagation failed
115     for (int ir=0;ir<fNBins;ir++) for (int ip=kNParTypes;ip--;) fData[ip][ir] += tmpAcc[ir*kNParTypes+ip];
116   }
117   //
118   for (int ip=kNParTypes;ip--;) for (int ir=fNBins;ir--;) fData[ip][ir] /= ntest;
119   delete[] tmpAcc;
120   //  Print();
121   //
122 }
123
124 //___________________________________________________________
125 void AliITSUMatLUT::Print(Option_t*) const
126 {
127   // print data 
128   printf("Average material budget in %d bins for %.4f<R<%.4f\n",fNBins,fRMin,fRMax);
129   printf("  # :  rMin : rMax  \t  X2X0  (  incr  )\t  RhoL  (  incr  )\n");
130   for (int i=0;i<fNBins;i++) {
131     double r = fRMin + fDR*i;
132     printf("%4d:%7.3f:%7.3f\t%.6f(%.6f)\t%.6f(%.6f)\n",i,r,r+fDR,
133            fData[kParX2X0][i],fData[kParX2X0][i]-(i==0 ? 0:fData[kParX2X0][i-1]),
134            fData[kParRhoL][i],fData[kParRhoL][i]-(i==0 ? 0:fData[kParRhoL][i-1]));
135   }
136
137 }
138
139 //___________________________________________________________
140 TH1* AliITSUMatLUT::GetHisto(const Option_t* option, const Char_t *name) const
141 {
142   // extract data to histo
143   if (fNBins<1) return 0;
144   TString nms = name;
145   if (nms.IsNull()) nms = "matLUT";
146   TH1F* h = new TH1F(nms.Data(),nms.Data(),fNBins,fRMin,fRMax);
147   TString opts = option;
148   opts.ToLower();
149   Bool_t diff = opts.Contains("d");
150   Int_t par = opts.Contains("rhol") ? kParRhoL : kParX2X0;
151   double valPrev = 0;
152   for (int i=0;i<fNBins;i++) {
153     double val = fData[par][i] - valPrev;
154     if (diff) valPrev = fData[par][i];
155     h->SetBinContent(i+1, val);
156   }
157   return h;
158 }
159
160 //___________________________________________________________
161 Double_t AliITSUMatLUT::GetMatBudget(const Double_t *pnt0, const Double_t *pnt1, Double_t *ret) const
162 {
163   // query the mat.budget between 2 points
164   double r0 = TMath::Sqrt(pnt0[0]*pnt0[0] + pnt0[1]*pnt0[1]);
165   double r1 = TMath::Sqrt(pnt1[0]*pnt1[0] + pnt1[1]*pnt1[1]);
166   if (r1<r0) {double t=r1;r1=r0;r0=t;}
167   double dr = r1-r0;
168   if (dr<1e-4) {
169     for (int i=kNParTypes;i--;) ret[i] = 0;
170     return 0;
171   }
172   //
173   double di = pnt1[0]-pnt0[0];
174   double dst = di*di;
175   di = pnt1[1]-pnt0[1];
176   dst += di*di;
177   di = pnt1[2]-pnt0[2];
178   dst += di*di;
179   dst = TMath::Sqrt(dst);
180   double dstf = dst/dr;
181   //
182   double par0[kNParTypes],par1[kNParTypes];
183   GetData(r0,par0);
184   GetData(r1,par1);
185   for (int i=kNParTypes;i--;) ret[i] = dstf*(par1[i]-par0[i]);
186   return dst;
187   //
188 }
189
190 //___________________________________________________________
191 void AliITSUMatLUT::GetData(Double_t r, Double_t* params) const
192 {
193   // get integrated parameters from rMin to r
194   double binF = (r-fRMin)*fDRInv;
195   if (binF<0) binF=0;
196   else if (binF>=fNBins) binF = fNBins-1;
197   int   bin   = int(binF);
198   double frac = 1.-(binF-bin);
199   for (int i=kNParTypes;i--;) {
200     double prev = bin ? fData[i][bin-1] : 0;
201     params[i] = fData[i][bin] - frac*(fData[i][bin]-prev);
202   }
203   //
204 }