ITS/UPGRADE/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   for (int itst=ntest;itst--;) {
  88     double parInt[kNParTypes]={0};
  89     double r   = fRMin;
  90     double phi = gRandom->Rndm()*TMath::Pi()*2;
  91     double cs  = TMath::Cos(phi);
  92     double sn  = TMath::Sin(phi);
  93     stop[0] = r*cs;
  94     stop[1] = r*sn;
  95     stop[2] = zmin + gRandom->Rndm()*(zmax-zmin);
  96     for (int ir=0;ir<fNBins;ir++) {
  97       r += dr;
  98       for (int i=3;i--;) start[i] = stop[i];
  99       stop[0] = r*cs;
 100       stop[1] = r*sn;
 101       AliTrackerBase::MeanMaterialBudget(start,stop, parStep);
 102       //
 103       parInt[kParX2X0] += parStep[1];
 104       parInt[kParRhoL] += parStep[0]*parStep[4];
 105       //
 106       for (int ip=kNParTypes;ip--;) fData[ip][ir] += parInt[ip];
 107     }
 108   }
 109   //
 110   for (int ip=kNParTypes;ip--;) for (int ir=fNBins;ir--;) fData[ip][ir] /= ntest;
 111   //
 112 }
 113
 114 //___________________________________________________________
 115 void AliITSUMatLUT::Print(Option_t*) const
 116 {
 117   // print data
 118   printf("Average material budget in %d bins for %.4f<R<%.4f\n",fNBins,fRMin,fRMax);
 119   printf("  # :  rMin : rMax  \t  X2X0  (  incr  )\t  RhoL  (  incr  )\n");
 120   for (int i=0;i<fNBins;i++) {
 121     double r = fRMin + fDR*i;
 122     printf("%4d:%7.3f:%7.3f\t%.6f(%.6f)\t%.6f(%.6f)\n",i,r,r+fDR,
 123            fData[kParX2X0][i],fData[kParX2X0][i]-(i==0 ? 0:fData[kParX2X0][i-1]),
 124            fData[kParRhoL][i],fData[kParRhoL][i]-(i==0 ? 0:fData[kParRhoL][i-1]));
 125   }
 126
 127 }
 128
 129 //___________________________________________________________
 130 TH1* AliITSUMatLUT::GetHisto(const Option_t* option, const Char_t *name) const
 131 {
 132   // extract data to histo
 133   if (fNBins<1) return 0;
 134   TString nms = name;
 135   if (nms.IsNull()) nms = "matLUT";
 136   TH1F* h = new TH1F(nms.Data(),nms.Data(),fNBins,fRMin,fRMax);
 137   TString opts = option;
 138   opts.ToLower();
 139   Bool_t diff = opts.Contains("d");
 140   Int_t par = opts.Contains("rhol") ? kParRhoL : kParX2X0;
 141   double valPrev = 0;
 142   for (int i=0;i<fNBins;i++) {
 143     double val = fData[par][i] - valPrev;
 144     if (diff) valPrev = fData[par][i];
 145     h->SetBinContent(i+1, val);
 146   }
 147   return h;
 148 }
 149
 150 //___________________________________________________________
 151 Double_t AliITSUMatLUT::GetMatBudget(const Double_t *pnt0, const Double_t *pnt1, Double_t *ret) const
 152 {
 153   // query the mat.budget between 2 points
 154   double r0 = TMath::Sqrt(pnt0[0]*pnt0[0] + pnt0[1]*pnt0[1]);
 155   double r1 = TMath::Sqrt(pnt1[0]*pnt1[0] + pnt1[1]*pnt1[1]);
 156   if (r1<r0) {double t=r1;r1=r0;r0=t;}
 157   double dr = r1-r0;
 158   if (dr<1e-4) {
 159     for (int i=kNParTypes;i--;) ret[i] = 0;
 160     return 0;
 161   }
 162   //
 163   double di = pnt1[0]-pnt0[0];
 164   double dst = di*di;
 165   di = pnt1[1]-pnt0[1];
 166   dst += di*di;
 167   di = pnt1[2]-pnt0[2];
 168   dst += di*di;
 169   dst = TMath::Sqrt(dst);
 170   double dstf = dst/dr;
 171   //
 172   double par0[kNParTypes],par1[kNParTypes];
 173   GetData(r0,par0);
 174   GetData(r1,par1);
 175   for (int i=kNParTypes;i--;) ret[i] = dstf*(par1[i]-par0[i]);
 176   return dst;
 177   //
 178 }
 179
 180 //___________________________________________________________
 181 void AliITSUMatLUT::GetData(Double_t r, Double_t* params) const
 182 {
 183   // get integrated parameters from rMin to r
 184   double binF = (r-fRMin)*fDRInv;
 185   if (binF<0) binF=0;
 186   else if (binF>=fNBins) binF = fNBins-1;
 187   int   bin   = int(binF);
 188   double frac = 1.-(binF-bin);
 189   for (int i=kNParTypes;i--;) {
 190     double prev = bin ? fData[i][bin-1] : 0;
 191     params[i] = fData[i][bin] - frac*(fData[i][bin]-prev);
 192   }
 193   //
 194 }