Added mat.budget estimation macro

diff --git a/ITS/UPGRADE/macros/matBudget.C b/ITS/UPGRADE/macros/matBudget.C
new file mode 100644 (file)
index 0000000..3825d52
--- /dev/null
+++ b/ITS/UPGRADE/macros/matBudget.C
@@ -0,0 +1,103 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "../ITS/UPGRADE/AliITSUTrackerGlo.h"
+#include "../ITS/UPGRADE/AliITSURecoDet.h"
+#include "../ITS/UPGRADE/AliITSURecoLayer.h"
+#include "TProfile.h"
+#include "TH2.h"
+#include "TAxis.h"
+#include "TRandom.h"
+#include "TFile.h"
+#endif
+
+/*
+  to run: 
+  .x LoadLibs.C 
+  .x chkTracker.C 
+  .x matBudget.C(trk)
+*/
+
+TProfile* histo[20]={0};
+TH2F *rphi=0, *rphic=0;
+
+void matBudget(AliITSUTrackerGlo* tracker, int nTest=100000, double phiMin=0,double phiMax=1.5708, int nbins=1000)
+{
+  AliITSURecoDet* its = tracker->GetITSInterface();
+  int nlr = its->GetNLayers();
+  double pnt0[3],pnt1[3],par[10];
+
+  for (int ilr=0;ilr<nlr;ilr++) {
+    histo[ilr] = new TProfile(Form("lr%d",ilr),
+                             Form("lr%d act%d",ilr,its->GetLayer(ilr)->GetActiveID()),
+                             nbins,phiMin,phiMax
+                             );
+  }
+  //
+  rphi  = new TH2F("rphi" ,"",nbins,phiMin,phiMax,its->GetRMax()*10,0,its->GetRMax());
+  rphic = new TH2F("rphic","",nbins,phiMin,phiMax,its->GetRMax()*10,0,its->GetRMax());
+  //
+  float frac = 0.1;
+  int npr=nTest*frac;
+  for (int it=0;it<nTest;it++) {
+    if ( (it%npr)==0 ) printf("Done %.1f%%\n",100*float(it)/nTest);
+    double phi = phiMin + gRandom->Rndm()*(phiMax-phiMin);
+    pnt1[0]=pnt1[1]=0;
+    pnt1[2]=gRandom->Rndm()*10;
+    double rmin=0;
+    for (int ilr=0;ilr<nlr;ilr++) {
+      double rmax = ilr<nlr-1 ? (its->GetLayer(ilr+1)->GetRMin()+its->GetLayer(ilr)->GetRMax())/2. : its->GetLayer(ilr)->GetRMax();
+      if (it==0) {
+       printf("Test of Lr%d in %f < R < %f\n",ilr,rmin,rmax);
+       histo[ilr]->SetTitle(Form("%s %.3f<R<%.3f",histo[ilr]->GetTitle(),rmin,rmax));
+      }
+      for (int i=3;i--;) pnt0[i]=pnt1[i];
+      pnt1[0] = rmax*TMath::Cos(phi);
+      pnt1[1] = rmax*TMath::Sin(phi);
+      pnt1[2] += rmax*1e-3; //to avoid strictly normal tracks
+      tracker->MeanMaterialBudget(pnt0,pnt1,par);
+      histo[ilr]->Fill(phi,par[1]);
+      rmin = rmax;
+    }   
+  }
+  //
+  for (int ilr=0;ilr<nlr;ilr++) {
+    if (!histo[ilr]->GetEntries()) continue;
+    histo[ilr]->Scale(100);
+    histo[ilr]->Fit("pol0","l");
+    double rmsy = histo[ilr]->GetRMS(2);
+    if (rmsy<1e-6) {
+      double cst = histo[ilr]->GetMean(2)*100;
+      histo[ilr]->SetMinimum(cst*0.8);      
+      histo[ilr]->SetMaximum(cst*1.2);
+    }
+  }
+  //
+  TAxis* xax = rphi->GetXaxis();
+  TAxis* yax = rphi->GetYaxis();
+  int nbr = yax->GetNbins();
+  for (int iphi=1;iphi<=nbins;iphi++) {
+    double phi = xax->GetBinCenter(iphi);
+    double rmax = 0;
+    pnt1[0]=pnt1[1]=0;
+    double avCum = 0;
+    for (int ir=2;ir<nbr;ir++) {
+      rmax = yax->GetBinCenter(ir);
+      pnt1[0] = rmax*TMath::Cos(phi);
+      pnt1[1] = rmax*TMath::Sin(phi);
+      for (int i=2;i--;) pnt0[i]=pnt1[i];
+      double av = 0;
+      for (int it=0;it<100;it++) {
+       pnt1[2] = pnt0[2] = gRandom->Rndm()*10;
+       pnt1[2] += rmax*1e-3; //to avoid strictly normal tracks
+       tracker->MeanMaterialBudget(pnt0,pnt1,par);
+       av += par[1];
+      }
+      av /= 10;
+      avCum += av;
+      rphi->SetBinContent(iphi,ir,av);
+      rphic->SetBinContent(iphi,ir,avCum);
+    }
+  }
+  //
+
+}