PWGUD/multVScentPbPb/ppcor.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include "TList.h"
   3 #include "TFile.h"
   4 #include "TStyle.h"
   5 #include "TH1F.h"
   6 #include "TH2F.h"
   7 #include "THnSparse.h"
   8 #include "TLegend.h"
   9 #include "TSystem.h"
  10 #include "TMath.h"
  11 #include "TCanvas.h"
  12 #include "TLegend.h"
  13 #include "TLatex.h"
  14 #include "TF1.h"
  15 #include "TLine.h"
  16 #include "TPaletteAxis.h"
  17 #include "TArrayD.h"
  18 #include "TGraphErrors.h"
  19 //
  20 //
  21 #endif
  22
  23
  24 TH2 *hMCAccTr,*hMCAccClsS,*hMCAccCls,*hMCGen,*hzvMC2;
  25 TH2 *hDtAccTr,*hDtAccClsS,*hDtAccCls,*hzvDt2;
  26 TH1 *hzvMC,*hzvDt,*hMCTrue;
  27 //
  28 TH2 *hTrAcceptance,*hTrDtAccNrm,*hTrDtCorr2D;
  29 TH2 *hClAcceptance,*hClDtAccNrm,*hClDtCorr2D;
  30 TH1 *hTrDtCorr1D,*hClDtCorr1D;
  31 //
  32 // for test of MC on MC only
  33 TH2 *hMCGenData = 0;
  34 TH1 *hMCTrueData = 0;
  35
  36 TH1* GetHisto(TList* lst, const char* name, const char* pref="");
  37 void CorrectMult(TH2* hMCAcc, TH2* hMCGn,
  38                  TH2* hDtAcc, TH1* hzDt,
  39                  const char* pref,
  40                  // output
  41                  TH2 *&hDtAccNrm,    // raw data normalized per Z bin
  42                  TH2 *&hAcceptance,  // 2D acceptance
  43                  TH2 *&hDtCorr2D,    // 2D corrected dN/dEta per Z bin (norm up to bin width)
  44                  TH1 *&hDtCorr1D     // 1D corrected dN/dEta
  45                  );
  46
  47
  48 const double minZStat = 50;
  49
  50 void ppcor(const char* flMC, const char* flDt)
  51 {
  52   TFile *fileMC = TFile::Open(flMC);
  53   TList* clistMC = (TList*)fileMC->Get("clist");
  54   TFile *fileDt = TFile::Open(flDt);
  55   TList* clistDt = (TList*)fileDt->Get("clist");
  56   //
  57   hMCAccTr   = (TH2*)GetHisto(clistMC,"b0_TrData_ZvEtaCutT","_mc");
  58   hMCAccClsS = (TH2*)GetHisto(clistMC,"b0_TrData_ZvEtaSPD1","_mc");
  59   hMCAccCls  = (TH2*)hMCAccClsS->Clone( Form("MCClusSPD_%s","_mc") );
  60   hMCAccCls->Add(hMCAccTr);
  61   hMCGen   = (TH2*)GetHisto(clistMC,"b0_zvEtaPrimMC","_mc");
  62   hzvMC2   = (TH2*)GetHisto(clistMC,"zv","_mc");
  63   hzvMC    =  hzvMC2->ProjectionX("zvgenMC");
  64   //
  65   hMCTrue = hMCGen->ProjectionX("MCTrue");
  66   if (hzvMC2->Integral()>0) hMCTrue->Scale(1./hzvMC2->Integral());
  67   hMCTrue->Scale(1./hMCTrue->GetBinWidth(1));
  68   //
  69   hDtAccTr   = (TH2*)GetHisto(clistDt,"b0_TrData_ZvEtaCutT","_dt");
  70   hDtAccClsS = (TH2*)GetHisto(clistDt,"b0_TrData_ZvEtaSPD1","_dt");
  71   hDtAccCls  = (TH2*)hDtAccClsS->Clone( Form("DtClusSPD_%s","_dt") );
  72   hDtAccCls->Add(hDtAccTr);
  73   hzvDt2   = (TH2*)GetHisto(clistDt,"zv","_dt");
  74   hzvDt   =  hzvDt2->ProjectionX("zvDt");
  75   //
  76   hMCGenData = (TH2*)GetHisto(clistDt,"b0_zvEtaPrimMC","_data");
  77   if (hMCGenData) {
  78     hMCTrueData = hMCGenData->ProjectionX("MCTrueData");
  79     if (hzvDt2->Integral()>0) hMCTrueData->Scale(1./hzvDt2->Integral());
  80     hMCTrueData->Scale(1./hMCTrueData->GetBinWidth(1));
  81   }
  82   //
  83   CorrectMult(hMCAccTr,hMCGen,
  84               hDtAccTr,hzvDt, "trc",
  85               hTrDtAccNrm,hTrAcceptance,hTrDtCorr2D,hTrDtCorr1D);
  86   //
  87   CorrectMult(hMCAccCls,hMCGen,
  88               hDtAccCls,hzvDt, "cls",
  89               hClDtAccNrm,hClAcceptance,hClDtCorr2D,hClDtCorr1D);
  90   //
  91   //
  92   hTrDtCorr1D->Draw();  // dndeta from tracklets
  93   hClDtCorr1D->Draw();  // dndeta from clusters
  94 }
  95
  96
  97 TH1* GetHisto(TList* lst, const char* name, const char* pref)
  98 {
  99   TH1* hst = (TH1*)lst->FindObject(name);
 100   if (!hst) return 0;
 101   TH1* hcl = (TH1*) hst->Clone( Form("%s%ss",hst->GetName(),pref) );
 102   return hcl;
 103 }
 104
 105 void CorrectMult(TH2* hMCAcc, TH2* hMCGn,
 106                  TH2* hDtAcc, TH1* hzDt,
 107                  const char* pref,
 108                  // output
 109                  TH2 *&hDtAccNrm,    // raw data normalized per Z bin
 110                  TH2 *&hAcceptance,  // 2D acceptance
 111                  TH2 *&hDtCorr2D,    // 2D corrected dN/dEta per Z bin (norm up to bin width)
 112                  TH1 *&hDtCorr1D     // 1D corrected dN/dEta
 113                  )
 114 {
 115   //
 116   int neta = hDtAcc->GetXaxis()->GetNbins();
 117   int nz   = hDtAcc->GetYaxis()->GetNbins();
 118   //
 119   hDtAccNrm = (TH2*) hDtAcc->Clone( Form("rawDataNormZ_%s",pref) );
 120   hDtAccNrm->Reset();
 121   for (int iz=1;iz<=nz;iz++) { // normalize data histo per Z bin
 122     double zv = hDtAcc->GetYaxis()->GetBinCenter(iz);
 123     int iz0 = hzDt->FindBin(zv);
 124     double zst  = hzDt->GetBinContent(iz0);
 125     double zstE = hzDt->GetBinError(iz0);
 126     if (zst<minZStat) continue;
 127     for (int ix=1;ix<=neta;ix++) {
 128       double vl = hDtAcc->GetBinContent(ix,iz);
 129       double er = hDtAcc->GetBinError(ix,iz);
 130       if (vl<1e-9 || er<1e-6) continue;
 131       double rat  = vl/zst;
 132       double ratE = rat*TMath::Sqrt( er*er/(vl*vl) + (zstE*zstE)/(zst*zst) );
 133       //printf("%2d %2d  %+e(%+e) -> %+e(%+e)\n",ix,iz,vl,er,rat,ratE);
 134       hDtAccNrm->SetBinContent(ix,iz,rat);
 135       hDtAccNrm->SetBinError(ix,iz,ratE);
 136     }
 137   }
 138   //
 139   hAcceptance = (TH2*)hMCAcc->Clone( Form("Acceptance_%s",pref) );
 140   hAcceptance->Divide(hMCGn);
 141   //
 142   hDtCorr2D = (TH2*)hDtAccNrm->Clone( Form("dNdEtaCor2D_%s",pref) );
 143   hDtCorr2D->Divide(hAcceptance);
 144   //
 145   hDtCorr1D = hDtCorr2D->ProjectionX( Form("dNdEtaCor1D_%s",pref) );
 146   hDtCorr1D->Reset();
 147   //
 148   for (int ix=1;ix<=neta;ix++) { // get weighted average
 149     double sm=0,sme=0;
 150     double bw = hDtCorr1D->GetBinWidth(ix);
 151     for (int iz=1;iz<=nz;iz++) {
 152       double vl = hDtCorr2D->GetBinContent(ix,iz);
 153       double er = hDtCorr2D->GetBinError(ix,iz);
 154       if (vl<1e-6 || er<1e-12) continue;
 155       sm += vl/(er*er);
 156       sme += 1./(er*er);
 157     }
 158     if (sme<1e-6) continue;
 159     sm /= sme;
 160     sme = 1./TMath::Sqrt(sme);
 161     hDtCorr1D->SetBinContent(ix,sm/bw);
 162     hDtCorr1D->SetBinError(ix,sme/bw);
 163   }
 164   //
 165 }