ROOT logo
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "TList.h"
#include "TFile.h"
#include "TStyle.h"
#include "TH1F.h"
#include "TH2F.h"
#include "THnSparse.h"
#include "TLegend.h"
#include "TSystem.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TF1.h"
#include "TLine.h"
#include "TPaletteAxis.h"
#include "TArrayD.h"
#include "TGraphErrors.h"
//
//
#endif


TH2 *hMCAccTr,*hMCAccClsS,*hMCAccCls,*hMCGen,*hzvMC2;
TH2 *hDtAccTr,*hDtAccClsS,*hDtAccCls,*hzvDt2;
TH1 *hzvMC,*hzvDt,*hMCTrue;
//
TH2 *hTrAcceptance,*hTrDtAccNrm,*hTrDtCorr2D;
TH2 *hClAcceptance,*hClDtAccNrm,*hClDtCorr2D;
TH1 *hTrDtCorr1D,*hClDtCorr1D;
//
// for test of MC on MC only
TH2 *hMCGenData = 0;
TH1 *hMCTrueData = 0;

TH1* GetHisto(TList* lst, const char* name, const char* pref="");
void CorrectMult(TH2* hMCAcc, TH2* hMCGn,
		 TH2* hDtAcc, TH1* hzDt,
		 const char* pref,
		 // output
		 TH2 *&hDtAccNrm,    // raw data normalized per Z bin
		 TH2 *&hAcceptance,  // 2D acceptance
		 TH2 *&hDtCorr2D,    // 2D corrected dN/dEta per Z bin (norm up to bin width)
		 TH1 *&hDtCorr1D     // 1D corrected dN/dEta 
		 );


const double minZStat = 50;

void ppcor(const char* flMC, const char* flDt)
{
  TFile *fileMC = TFile::Open(flMC);
  TList* clistMC = (TList*)fileMC->Get("clist");
  TFile *fileDt = TFile::Open(flDt);
  TList* clistDt = (TList*)fileDt->Get("clist");
  //
  hMCAccTr   = (TH2*)GetHisto(clistMC,"b0_TrData_ZvEtaCutT","_mc"); 
  hMCAccClsS = (TH2*)GetHisto(clistMC,"b0_TrData_ZvEtaSPD1","_mc");
  hMCAccCls  = (TH2*)hMCAccClsS->Clone( Form("MCClusSPD_%s","_mc") );
  hMCAccCls->Add(hMCAccTr);
  hMCGen   = (TH2*)GetHisto(clistMC,"b0_zvEtaPrimMC","_mc");
  hzvMC2   = (TH2*)GetHisto(clistMC,"zv","_mc");
  hzvMC    =  hzvMC2->ProjectionX("zvgenMC");
  //
  hMCTrue = hMCGen->ProjectionX("MCTrue");
  if (hzvMC2->Integral()>0) hMCTrue->Scale(1./hzvMC2->Integral());
  hMCTrue->Scale(1./hMCTrue->GetBinWidth(1));
  //
  hDtAccTr   = (TH2*)GetHisto(clistDt,"b0_TrData_ZvEtaCutT","_dt"); 
  hDtAccClsS = (TH2*)GetHisto(clistDt,"b0_TrData_ZvEtaSPD1","_dt"); 
  hDtAccCls  = (TH2*)hDtAccClsS->Clone( Form("DtClusSPD_%s","_dt") );
  hDtAccCls->Add(hDtAccTr);
  hzvDt2   = (TH2*)GetHisto(clistDt,"zv","_dt");
  hzvDt   =  hzvDt2->ProjectionX("zvDt");
  //
  hMCGenData = (TH2*)GetHisto(clistDt,"b0_zvEtaPrimMC","_data");
  if (hMCGenData) {
    hMCTrueData = hMCGenData->ProjectionX("MCTrueData");
    if (hzvDt2->Integral()>0) hMCTrueData->Scale(1./hzvDt2->Integral());
    hMCTrueData->Scale(1./hMCTrueData->GetBinWidth(1));
  }
  //
  CorrectMult(hMCAccTr,hMCGen,
	      hDtAccTr,hzvDt, "trc",
	      hTrDtAccNrm,hTrAcceptance,hTrDtCorr2D,hTrDtCorr1D);
  //
  CorrectMult(hMCAccCls,hMCGen,
	      hDtAccCls,hzvDt, "cls",
	      hClDtAccNrm,hClAcceptance,hClDtCorr2D,hClDtCorr1D);
  //
  //
  hTrDtCorr1D->Draw();  // dndeta from tracklets
  hClDtCorr1D->Draw();  // dndeta from clusters
}


TH1* GetHisto(TList* lst, const char* name, const char* pref)
{
  TH1* hst = (TH1*)lst->FindObject(name);
  if (!hst) return 0;
  TH1* hcl = (TH1*) hst->Clone( Form("%s%ss",hst->GetName(),pref) );
  return hcl;
}

void CorrectMult(TH2* hMCAcc, TH2* hMCGn,
		 TH2* hDtAcc, TH1* hzDt,
		 const char* pref,
		 // output
		 TH2 *&hDtAccNrm,    // raw data normalized per Z bin
		 TH2 *&hAcceptance,  // 2D acceptance
		 TH2 *&hDtCorr2D,    // 2D corrected dN/dEta per Z bin (norm up to bin width)
		 TH1 *&hDtCorr1D     // 1D corrected dN/dEta 
		 )
{
  //
  int neta = hDtAcc->GetXaxis()->GetNbins();
  int nz   = hDtAcc->GetYaxis()->GetNbins();  
  //
  hDtAccNrm = (TH2*) hDtAcc->Clone( Form("rawDataNormZ_%s",pref) );
  hDtAccNrm->Reset();
  for (int iz=1;iz<=nz;iz++) { // normalize data histo per Z bin
    double zv = hDtAcc->GetYaxis()->GetBinCenter(iz);
    int iz0 = hzDt->FindBin(zv);
    double zst  = hzDt->GetBinContent(iz0);
    double zstE = hzDt->GetBinError(iz0);
    if (zst<minZStat) continue;
    for (int ix=1;ix<=neta;ix++) {
      double vl = hDtAcc->GetBinContent(ix,iz);
      double er = hDtAcc->GetBinError(ix,iz);
      if (vl<1e-9 || er<1e-6) continue;
      double rat  = vl/zst;
      double ratE = rat*TMath::Sqrt( er*er/(vl*vl) + (zstE*zstE)/(zst*zst) );
      //printf("%2d %2d  %+e(%+e) -> %+e(%+e)\n",ix,iz,vl,er,rat,ratE);
      hDtAccNrm->SetBinContent(ix,iz,rat);
      hDtAccNrm->SetBinError(ix,iz,ratE);
    }
  }
  //
  hAcceptance = (TH2*)hMCAcc->Clone( Form("Acceptance_%s",pref) );
  hAcceptance->Divide(hMCGn);
  //
  hDtCorr2D = (TH2*)hDtAccNrm->Clone( Form("dNdEtaCor2D_%s",pref) );
  hDtCorr2D->Divide(hAcceptance);
  //
  hDtCorr1D = hDtCorr2D->ProjectionX( Form("dNdEtaCor1D_%s",pref) );
  hDtCorr1D->Reset();
  //
  for (int ix=1;ix<=neta;ix++) { // get weighted average
    double sm=0,sme=0;
    double bw = hDtCorr1D->GetBinWidth(ix);
    for (int iz=1;iz<=nz;iz++) {
      double vl = hDtCorr2D->GetBinContent(ix,iz);
      double er = hDtCorr2D->GetBinError(ix,iz);
      if (vl<1e-6 || er<1e-12) continue;
      sm += vl/(er*er);
      sme += 1./(er*er);
    }
    if (sme<1e-6) continue;
    sm /= sme;
    sme = 1./TMath::Sqrt(sme);
    hDtCorr1D->SetBinContent(ix,sm/bw);
    hDtCorr1D->SetBinError(ix,sme/bw);
  }
  //
}
 ppcor.C:1
 ppcor.C:2
 ppcor.C:3
 ppcor.C:4
 ppcor.C:5
 ppcor.C:6
 ppcor.C:7
 ppcor.C:8
 ppcor.C:9
 ppcor.C:10
 ppcor.C:11
 ppcor.C:12
 ppcor.C:13
 ppcor.C:14
 ppcor.C:15
 ppcor.C:16
 ppcor.C:17
 ppcor.C:18
 ppcor.C:19
 ppcor.C:20
 ppcor.C:21
 ppcor.C:22
 ppcor.C:23
 ppcor.C:24
 ppcor.C:25
 ppcor.C:26
 ppcor.C:27
 ppcor.C:28
 ppcor.C:29
 ppcor.C:30
 ppcor.C:31
 ppcor.C:32
 ppcor.C:33
 ppcor.C:34
 ppcor.C:35
 ppcor.C:36
 ppcor.C:37
 ppcor.C:38
 ppcor.C:39
 ppcor.C:40
 ppcor.C:41
 ppcor.C:42
 ppcor.C:43
 ppcor.C:44
 ppcor.C:45
 ppcor.C:46
 ppcor.C:47
 ppcor.C:48
 ppcor.C:49
 ppcor.C:50
 ppcor.C:51
 ppcor.C:52
 ppcor.C:53
 ppcor.C:54
 ppcor.C:55
 ppcor.C:56
 ppcor.C:57
 ppcor.C:58
 ppcor.C:59
 ppcor.C:60
 ppcor.C:61
 ppcor.C:62
 ppcor.C:63
 ppcor.C:64
 ppcor.C:65
 ppcor.C:66
 ppcor.C:67
 ppcor.C:68
 ppcor.C:69
 ppcor.C:70
 ppcor.C:71
 ppcor.C:72
 ppcor.C:73
 ppcor.C:74
 ppcor.C:75
 ppcor.C:76
 ppcor.C:77
 ppcor.C:78
 ppcor.C:79
 ppcor.C:80
 ppcor.C:81
 ppcor.C:82
 ppcor.C:83
 ppcor.C:84
 ppcor.C:85
 ppcor.C:86
 ppcor.C:87
 ppcor.C:88
 ppcor.C:89
 ppcor.C:90
 ppcor.C:91
 ppcor.C:92
 ppcor.C:93
 ppcor.C:94
 ppcor.C:95
 ppcor.C:96
 ppcor.C:97
 ppcor.C:98
 ppcor.C:99
 ppcor.C:100
 ppcor.C:101
 ppcor.C:102
 ppcor.C:103
 ppcor.C:104
 ppcor.C:105
 ppcor.C:106
 ppcor.C:107
 ppcor.C:108
 ppcor.C:109
 ppcor.C:110
 ppcor.C:111
 ppcor.C:112
 ppcor.C:113
 ppcor.C:114
 ppcor.C:115
 ppcor.C:116
 ppcor.C:117
 ppcor.C:118
 ppcor.C:119
 ppcor.C:120
 ppcor.C:121
 ppcor.C:122
 ppcor.C:123
 ppcor.C:124
 ppcor.C:125
 ppcor.C:126
 ppcor.C:127
 ppcor.C:128
 ppcor.C:129
 ppcor.C:130
 ppcor.C:131
 ppcor.C:132
 ppcor.C:133
 ppcor.C:134
 ppcor.C:135
 ppcor.C:136
 ppcor.C:137
 ppcor.C:138
 ppcor.C:139
 ppcor.C:140
 ppcor.C:141
 ppcor.C:142
 ppcor.C:143
 ppcor.C:144
 ppcor.C:145
 ppcor.C:146
 ppcor.C:147
 ppcor.C:148
 ppcor.C:149
 ppcor.C:150
 ppcor.C:151
 ppcor.C:152
 ppcor.C:153
 ppcor.C:154
 ppcor.C:155
 ppcor.C:156
 ppcor.C:157
 ppcor.C:158
 ppcor.C:159
 ppcor.C:160
 ppcor.C:161
 ppcor.C:162
 ppcor.C:163
 ppcor.C:164
 ppcor.C:165
 ppcor.C:166