1 #include "AliITSUMatLUT.h"
3 #include "AliTrackerBase.h"
10 //___________________________________________________________
11 AliITSUMatLUT::AliITSUMatLUT()
19 for (int i=kNParTypes;i--;) fData[i]=0;
22 //___________________________________________________________
23 AliITSUMatLUT::~AliITSUMatLUT()
26 for (int i=kNParTypes;i--;) delete fData[i];
29 //___________________________________________________________
30 AliITSUMatLUT::AliITSUMatLUT(Double_t rmin,Double_t rmax,Int_t nbin)
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));
48 //___________________________________________________________
49 AliITSUMatLUT::AliITSUMatLUT(const AliITSUMatLUT& src)
58 for (int i=kNParTypes;i--;) {
60 fData[i] = new Double_t[fNBins];
61 memcpy(fData[i],src.fData[i],fNBins*sizeof(Double_t));
67 //___________________________________________________________
68 AliITSUMatLUT & AliITSUMatLUT::operator=(const AliITSUMatLUT& src)
71 if (this == &src) return *this;
72 this->~AliITSUMatLUT();
73 new(this) AliITSUMatLUT(src);
78 //___________________________________________________________
79 void AliITSUMatLUT::FillData(Int_t ntest, Double_t zmin,Double_t zmax)
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};
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;
98 stop[2] = zmin + gRandom->Rndm()*(zmax-zmin);
100 for (int ir=0;ir<fNBins;ir++) {
102 for (int i=3;i--;) start[i] = stop[i];
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;}
109 parInt[kParX2X0] += parStep[1];
110 parInt[kParRhoL] += parStep[0]*parStep[4];
112 for (int ip=kNParTypes;ip--;) tmpAcc[ir*kNParTypes+ip] = parInt[ip];
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];
118 for (int ip=kNParTypes;ip--;) for (int ir=fNBins;ir--;) fData[ip][ir] /= ntest;
124 //___________________________________________________________
125 void AliITSUMatLUT::Print(Option_t*) const
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]));
139 //___________________________________________________________
140 TH1* AliITSUMatLUT::GetHisto(const Option_t* option, const Char_t *name) const
142 // extract data to histo
143 if (fNBins<1) return 0;
145 if (nms.IsNull()) nms = "matLUT";
146 TH1F* h = new TH1F(nms.Data(),nms.Data(),fNBins,fRMin,fRMax);
147 TString opts = option;
149 Bool_t diff = opts.Contains("d");
150 Int_t par = opts.Contains("rhol") ? kParRhoL : kParX2X0;
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);
160 //___________________________________________________________
161 Double_t AliITSUMatLUT::GetMatBudget(const Double_t *pnt0, const Double_t *pnt1, Double_t *ret) const
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;}
169 for (int i=kNParTypes;i--;) ret[i] = 0;
173 double di = pnt1[0]-pnt0[0];
175 di = pnt1[1]-pnt0[1];
177 di = pnt1[2]-pnt0[2];
179 dst = TMath::Sqrt(dst);
180 double dstf = dst/dr;
182 double par0[kNParTypes],par1[kNParTypes];
185 for (int i=kNParTypes;i--;) ret[i] = dstf*(par1[i]-par0[i]);
190 //___________________________________________________________
191 void AliITSUMatLUT::GetData(Double_t r, Double_t* params) const
193 // get integrated parameters from rMin to r
194 double binF = (r-fRMin)*fDRInv;
196 else if (binF>=fNBins) binF = fNBins-1;
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);