ROOT logo
/*************************************************************************
* Macro correctSPDdNdEta                                                 *
* To read data and correction histograms produced with                   *
* AliAnalysisTaskSPDdNdEta and apply corrections to data                 *    
*                                                                        *
* Author:  M. Nicassio (INFN Bari)                                       *
* Contact: Maria.Nicassio@ba.infn.it, Domenico.Elia@ba.infn.it           *
**************************************************************************/

#include "Riostream.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TList.h" 
#include "TTree.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TCanvas.h"

Float_t CalculateErrordNdEtaPerPartPair(Float_t a, Float_t b, Float_t da, Float_t db); 
void combscalingfactors (TList *lall, TList *lcomb, TList *lmc, Double_t *scaleNormRange_rot, Double_t *scaleNormRange_labels);

void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, 
                       Char_t* fileAlphaCorr, Char_t* fileMCBkgCorr, Char_t* fileMC, 
                       Int_t binVtxStart=5, Int_t binVtxStop=15,
                       Double_t cutCorr=10, Double_t cutBkg=1., Double_t cutBkgMC=1., 
                       Bool_t bkgFromMCLabels = kFALSE,
                       Bool_t changeStrangeness = kFALSE)   { 

// vtx lim 13-27 for +-7cm
// vtx lim 15-25 for +-5cm


  Int_t nBinVtx = 40;
  Double_t MaxVtx = 20.;

  const Int_t nBinMultCorr = 200; 
  Float_t nMaxMult = 20000.;
  Double_t binLimMultCorr[nBinMultCorr+1];
  binLimMultCorr[0] = 0.;
  binLimMultCorr[1] = 1.;
  for (Int_t i = 2; i<=nBinMultCorr;++i) {
    binLimMultCorr[i] = (i-1)*nMaxMult/nBinMultCorr;
  }

  const Int_t nBinEtaCorr = 60; 
  Float_t etaMax = 3.;
  Float_t etaMin = -3.;
  Double_t *binLimEtaCorr = new Double_t[nBinEtaCorr+1];
  for (Int_t i = 0; i<nBinEtaCorr+1; ++i) {
    binLimEtaCorr[i] = (Double_t) etaMin+i*(etaMax*2.)/nBinEtaCorr;
  }
  
  Float_t nBinsPerPseudorapidityUnit = nBinEtaCorr/(binLimEtaCorr[nBinEtaCorr]-binLimEtaCorr[0]);
  cout<<"Number of bins per pseudorapidity unit"<<nBinsPerPseudorapidityUnit<<endl;

  gStyle->SetOptLogy(kFALSE);
  gStyle->SetFrameLineWidth(2);
  gStyle->SetFrameFillColor(0); 
  gStyle->SetCanvasColor(0);
  gStyle->SetPadColor(0);
  gStyle->SetHistLineWidth(2);
  gStyle->SetLabelSize(0.05, "x");
  gStyle->SetLabelSize(0.05, "y");
  gStyle->SetTitleSize(0.05, "x");
  gStyle->SetTitleSize(0.05, "y");
  gStyle->SetTitleOffset(1.1, "x");
  gStyle->SetTitleOffset(0.9, "y"); 
  gStyle->SetPadBottomMargin(0.14);
  gStyle->SetPalette(1,0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetStatColor(0);

  //Input files
  TFile *f_MC = new TFile(fileMC);
  TFile *f_acorr = new TFile(fileAlphaCorr);
  TFile *f_mcbkgcorr = new TFile(fileMCBkgCorr);
  TFile *f_rawbkg = new TFile(fileRawBkgCorr);
  TFile *f_raw = new TFile(fileRaw);

  TList *l_MC = (TList*)f_MC->Get("cOutput");
  TList *l_acorr = (TList*)f_acorr->Get("cOutput");
  TList *l_mcbkgcorr = (TList*)f_mcbkgcorr->Get("cOutput");
  TList *l_rawbkg = (TList*)f_rawbkg->Get("cOutput");
  TList *l_raw = (TList*)f_raw->Get("cOutput");

  Double_t scaleBkg = 1.;
  Double_t scaleBkgMC = 1.;
  Double_t scaleBkgLabels = 1.;

  cout<<" Background scaling factors for MC: "<<endl;
  combscalingfactors (l_acorr,l_mcbkgcorr,l_acorr,&scaleBkgMC,&scaleBkgLabels);
  cout<<" Background scaling factors for data: "<<endl;
  combscalingfactors (l_raw,l_rawbkg,l_acorr,&scaleBkg,&scaleBkgLabels);

  //Histogram to be corrected at tracklet level
  TH2F *hSPDEta_2Draw = new TH2F(*((TH2F*)l_raw->FindObject("fHistSPDRAWEtavsZ")));
  hSPDEta_2Draw->SetNameTitle("SPDEta_2Draw","Reconstructed tracklets");

  //Corrections at track level
  // Combinatorial background: beta correction from data
  TH2F *hCombBkg = (TH2F*)l_rawbkg->FindObject("fHistSPDRAWEtavsZ");
  hCombBkg->Scale(scaleBkg);

  TH2F *hCombBkgCorrData = new TH2F("backgroundCorrDATA","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hCombBkgCorrData->GetXaxis()->SetTitle("#eta");
  hCombBkgCorrData->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
  hCombBkgCorrData->Sumw2();
  
  // Combinatorial background: beta correction from MC
  TH2F *hCombBkgMC = (TH2F*)l_mcbkgcorr->FindObject("fHistSPDRAWEtavsZ"); 
  hCombBkgMC->Scale(scaleBkgMC);
  TH2F *hCombBkgMCDen = (TH2F*)l_acorr->FindObject("fHistSPDRAWEtavsZ");

  TH2F *hCombBkgCorrMC = new TH2F("backgroundCorrMC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hCombBkgCorrMC->GetXaxis()->SetTitle("#eta");
  hCombBkgCorrMC->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
  hCombBkgCorrMC->Sumw2();
 
  // Background correction from MC labels
  TH2F *hBkgCombLabels = (TH2F*)l_acorr->FindObject("fHistBkgCombLabels");
  if (bkgFromMCLabels) {
    hCombBkgCorrData->Divide(hBkgCombLabels,hSPDEta_2Draw,1.,1.);
    hCombBkgCorrData->Scale(scaleBkgLabels);
  } else {
    hCombBkgCorrData->Divide(hCombBkg,hSPDEta_2Draw,1.,1.);
  }

  // 1-beta in DATA
  TH2F *hCombBkgCorr2Data = new TH2F("backgroundCorr2Data","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);

  for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
    for (Int_t jz=0; jz<nBinVtx; jz++) {
      hCombBkgCorr2Data->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrData->GetBinContent(ieta+1,jz+1));
      hCombBkgCorr2Data->SetBinError(ieta+1,jz+1,hCombBkgCorrData->GetBinError(ieta+1,jz+1));
    }
  }
 
  // Errors on beta
  //..data
  for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
    for (Int_t jz=0; jz<nBinVtx; jz++) {
      hCombBkgCorrData->SetBinError(ieta+1,jz+1,scaleBkg*TMath::Sqrt(hCombBkg->GetBinContent(ieta+1,jz+1))/hSPDEta_2Draw->GetBinContent(ieta+1,jz+1));
//      hCombBkgCorrData->SetBinError(ieta+1,jz+1,hCombBkgCorrData->GetBincontent(ieta+1,jz+1)*TMath::Sqrt(hSPDEta_2Draw->GetBinContent(ieta+1,jz+1))/hSPDEta_2Draw->GetBinContent(ieta+1,jz+1));
    }
  }

  // Alpha correction: MC only
  TH2F *hPrimaries_evtTrVtx = new TH2F(*((TH2F*) l_acorr->FindObject("fHistNonDetectableCorrNum"))); // all prim in sel events
  new TCanvas();
  hPrimaries_evtTrVtx->Draw(); 
 
  TH2F *hInvAlphaCorr = new TH2F("InvAlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hInvAlphaCorr->Sumw2();

  TH2F *hGenProtons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistProtons")));
  TH2F *hGenKaons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistKaons")));
  TH2F *hGenPions = new TH2F(*((TH2F*) l_acorr->FindObject("fHistPions")));
  TH2F *hGenOthers = new TH2F(*((TH2F*) l_acorr->FindObject("fHistOthers")));

  TH2F *hPrimReweighted = new TH2F("primReweighted","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hPrimReweighted->Add(hGenKaons);
  hPrimReweighted->Scale(2.25); //double kaon fraction
  hPrimReweighted->Add(hGenProtons);
  hPrimReweighted->Add(hGenPions);
  hPrimReweighted->Add(hGenOthers); //use this in alpha if flag for syst true

  TH2F *hRecProtons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedProtons")));
  TH2F *hRecKaons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedKaons")));
  TH2F *hRecPions = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedPions")));
  TH2F *hRecOthers = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedOthers")));
  TH2F *hRecSec = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedSec")));
 
  TH2F *hRecReweighted = new TH2F("recReweighted","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hRecReweighted->Add(hRecKaons);
  hRecReweighted->Scale(2.25);
  hRecReweighted->Add(hRecProtons);
  hRecReweighted->Add(hRecPions);
  hRecReweighted->Add(hRecOthers); //use this in alpha if flag for syst true
  hRecReweighted->Add(hRecSec);
  hRecReweighted->Add(hBkgCombLabels);

  // 1-beta in MC
  TH2F *hCombBkgCorr2MC = new TH2F("backgroundCorr2MC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);

  if (changeStrangeness) {
    if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hRecReweighted,1.,1.);
    else hCombBkgCorrMC->Divide(hCombBkgMC,hRecReweighted,1.,1.);

    for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
      for (Int_t jz=0; jz<nBinVtx; jz++)
        hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrMC->GetBinContent(ieta+1,jz+1));

    hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hRecReweighted);
    hInvAlphaCorr->Divide(hPrimReweighted);
  } else {

    if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hCombBkgMCDen,1.,1.);
    else hCombBkgCorrMC->Divide(hCombBkgMC,hCombBkgMCDen,1.,1.);

    for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
      for (Int_t jz=0; jz<nBinVtx; jz++)
        hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrMC->GetBinContent(ieta+1,jz+1));

    hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hCombBkgMCDen);
    hInvAlphaCorr->Divide(hPrimaries_evtTrVtx);

  }

//  new TCanvas();
//  hCombBkgCorr2MC->Draw();


  // Efficiency for particle species
  TH2F *hEffProtons = new TH2F("effProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  TH2F *hEffKaons   = new TH2F("effKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  TH2F *hEffPions   = new TH2F("effPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  TH2F *hEffOthers  = new TH2F("effOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);

  hEffProtons->Divide(hRecProtons,hGenProtons,1.,1.);  
  hEffKaons->Divide(hRecKaons,hGenKaons,1.,1.);
  hEffPions->Divide(hRecPions,hGenPions,1.,1.);
  hEffOthers->Divide(hRecOthers,hGenOthers,1.,1.);

//  hEffKaons->Divide(hEffPions);
  

  new TCanvas();
  hInvAlphaCorr->Draw();

  // Errors on alpha

  for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
    for (Int_t jz=0; jz<nBinVtx; jz++) {
      hInvAlphaCorr->SetBinError(ieta+1,jz+1,TMath::Sqrt((hCombBkgCorr2MC->GetBinContent(ieta+1,jz+1))*(hCombBkgMCDen->GetBinContent(ieta+1,jz+1)))/hPrimaries_evtTrVtx->GetBinContent(ieta+1,jz+1));
    }
  }

  TH2F *hAlphaCorr = new TH2F("AlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hAlphaCorr->GetXaxis()->SetTitle("#eta");
  hAlphaCorr->GetYaxis()->SetTitle("z_{vtx} [cm]");

  for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
    for (Int_t jz=0; jz<nBinVtx; jz++) {
       hAlphaCorr->SetBinContent(ieta+1,jz+1,1./hInvAlphaCorr->GetBinContent(ieta+1,jz+1));
    }
  }

  new TCanvas();
  hAlphaCorr->Draw();

 
  //////////////////////// Factorize alpha ///////////////////////////
   
  //////////////////////////////////////////////////////////////////// 

  // Number of events and normalization histograms
  TH1D *hSPDvertex = (TH1D*)l_raw->FindObject("fHistSPDvtxAnalysis");
  Double_t nEvtsRec=hSPDvertex->Integral(binVtxStart+1,binVtxStop);
  cout<<"Number of reconstructed events in the selected vertex range: "<<nEvtsRec <<endl;

  //Cutting corrections at the edges of the acceptance
  for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
    for (Int_t kvtx=0; kvtx<nBinVtx; kvtx++) {
      if (hAlphaCorr->GetBinContent(jeta+1,kvtx+1)>cutCorr||hAlphaCorr->GetBinContent(jeta+1,kvtx+1)<0) {
        hAlphaCorr->SetBinContent(jeta+1,kvtx+1,0.);
        hAlphaCorr->SetBinError(jeta+1,kvtx+1,0.);
        hInvAlphaCorr->SetBinContent(jeta+1,kvtx+1,0.);
        hInvAlphaCorr->SetBinError(jeta+1,kvtx+1,0.);
      }
      if (hCombBkgCorrData->GetBinContent(jeta+1,kvtx+1)>cutBkg) { 
        hCombBkgCorrData->SetBinContent(jeta+1,kvtx+1,0.);
        hCombBkgCorrData->SetBinError(jeta+1,kvtx+1,0.);
      }
      if (hCombBkgCorrMC->GetBinContent(jeta+1,kvtx+1)>cutBkgMC) {
        hCombBkgCorrMC->SetBinContent(jeta+1,kvtx+1,0.);
        hCombBkgCorrMC->SetBinError(jeta+1,kvtx+1,0.);
      }
    }
  }

  new TCanvas();
  hCombBkgCorrMC->Draw();
  new TCanvas();
  hCombBkgCorrData->Draw();

  // MC distribution 

  //... to compare corrected distribution to MC in selected events
  TH2F *hVertexMC2D = (TH2F*)l_MC->FindObject("fHistSelEvts"); 
  TH1D *hMCvertex = new TH1D(*(hVertexMC2D->ProjectionY("MCvertex",0,-1,"e")));  //MC vertex distrib
  new TCanvas();
  hMCvertex->Draw(); 

  TH2F *Eta = (TH2F*)l_MC->FindObject("fHistNonDetectableCorrNum"); 

  TH1D *hdNdEta = new TH1D("dNdEta","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);
  hdNdEta = Eta->ProjectionX("dNdEta",0,-1,"e"); //here already all events
  Double_t nEvtsTot=hMCvertex->Integral(0,hMCvertex->GetNbinsX()+1);

  Double_t nEvtsTotSel= hMCvertex->Integral(binVtxStart+1,binVtxStop);

  cout<<"Number of generated events: "<<nEvtsTot<<endl;
  cout<<"Number of generated events in the selected vertex range: "<<nEvtsTotSel <<endl;
  
  Float_t nprimcentraleta = Eta->ProjectionX("dNdEta",0,-1,"e")->Integral(26,35); // project all or only in the sel vertex range? 
  hdNdEta->Scale(nBinsPerPseudorapidityUnit/nEvtsTot);  // if so change number of events
  hdNdEta->GetXaxis()->SetTitle("#eta");
  hdNdEta->GetYaxis()->SetTitle("dN_{ch}/d#eta");
  for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta->SetBinError(i+1,0.);
  // dNdEta in +- 0.5
  cout<<"Monte Carlo dN/dEta in |eta|<0.5  (selected events for analysis and all events!): "<<nprimcentraleta/nEvtsTot<<endl;  
  cout<<"Monte Carlo dN/dEta in |eta|<0.5  (selected events for analysis and events in vtx range!): "<<(Eta->ProjectionX("_x",binVtxStart+1,binVtxStop,"e")->Integral(26,35))/nEvtsTotSel<<endl;
  new TCanvas();
  hdNdEta->Draw();


  //Corrected distributions
  TH2F *hSPDEta_2D_bkgCorr =  new TH2F("SPDEta_2D_bkgCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hSPDEta_2D_bkgCorr->Sumw2();

  TH2F *hSPDEta_2D_fullyCorr =  new TH2F("SPDEta_2D_fullyCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hSPDEta_2D_fullyCorr->Sumw2();

  TH1F *hnorm_fullyCorr =  new TH1F("norm_fullyCorr", "",nBinEtaCorr,binLimEtaCorr);

  //dN/deta
  TH1F *hdNdEta_raw = new TH1F("dNdEta_raw","",nBinEtaCorr,binLimEtaCorr);
  hdNdEta_raw->Sumw2();
  hdNdEta_raw->GetXaxis()->SetTitle("#eta");
  hdNdEta_raw->GetYaxis()->SetTitle("dN_{ch}/d#eta");
  hdNdEta_raw->Add(hSPDEta_2Draw->ProjectionX("SPDEta_2Draw_x",binVtxStart+1,binVtxStop,"e"));
  hdNdEta_raw->Scale(nBinsPerPseudorapidityUnit/nEvtsRec);
  new TCanvas();
  hdNdEta_raw->Draw(); 

  TH1F *hdNdEta_bkgCorr = new TH1F("dNdEta_bkgCorr","",nBinEtaCorr,binLimEtaCorr);
  hdNdEta_bkgCorr->Sumw2();
  TH1F *hdNdEta_fullyCorr = new TH1F("dNdEta_fullyCorr","",nBinEtaCorr,binLimEtaCorr);
  hdNdEta_fullyCorr->Sumw2();

  //Apply corrections at
  //...track level
  hSPDEta_2D_bkgCorr->Add(hSPDEta_2Draw); 
  hSPDEta_2D_bkgCorr->Multiply(hCombBkgCorr2Data);
  hSPDEta_2D_fullyCorr->Add(hSPDEta_2D_bkgCorr);
  hSPDEta_2D_fullyCorr->Divide(hInvAlphaCorr); 
  new TCanvas();
  hSPDEta_2D_bkgCorr->Draw();
  new TCanvas();
  hSPDEta_2D_fullyCorr->Draw();

  //Filling normalization histogram
  for (Int_t ivtx=binVtxStart; ivtx<binVtxStop; ivtx++) {
    Double_t nEvtsivtx = hSPDvertex->GetBinContent(ivtx+1);
    for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
      if (hAlphaCorr->GetBinContent(jeta+1,ivtx+1)!=0.) { 
        hnorm_fullyCorr->Fill(hAlphaCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtx);
      }
    }
  }
  
  for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
    hnorm_fullyCorr->SetBinError(jeta+1,0);
  }
  new TCanvas();
  hnorm_fullyCorr->Draw();

  hdNdEta_bkgCorr->Add(hSPDEta_2D_bkgCorr->ProjectionX("SPDEta_2D_bkgCorr_x",binVtxStart+1,binVtxStop,"e"));
  hdNdEta_bkgCorr->Scale(1./nEvtsRec); 
  hdNdEta_fullyCorr->Divide(hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e"),hnorm_fullyCorr);

  TH1F *hCorrecteddNdEtaCentr = new TH1F("correcteddNdEtaCentr","",10,-0.5,0.5);
  hCorrecteddNdEtaCentr->Sumw2();
  for (Int_t jeta=0; jeta<10; jeta++) {
    hCorrecteddNdEtaCentr->SetBinContent(jeta+1,hdNdEta_fullyCorr->GetBinContent(26+jeta));
    hCorrecteddNdEtaCentr->SetBinError(jeta+1,hdNdEta_fullyCorr->GetBinError(26+jeta));
  }
  hCorrecteddNdEtaCentr->Rebin(10);

  new TCanvas();
  hCorrecteddNdEtaCentr->Draw();

  Float_t ncorrtrackletscentr = hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e")->Integral(26,35);
  Float_t nevtscentr = hnorm_fullyCorr->Integral(26,35)/10.;
  cout<<""<<endl;
  cout<<"Corrected dN/dEta in |eta|<0.5: "<<ncorrtrackletscentr/nevtscentr<<endl;
  cout<<"Error on corrected tracklets in |eta|<0.5: +-"<<hCorrecteddNdEtaCentr->GetBinError(1)<<endl;
  // dN/dEta per participant pairs

  Float_t NpartV0Glauber = 381.188;
  Float_t NpartCl2Glauber = 378.5;
  Float_t NpartCl2Hij = 365.3;
  
  Float_t errNpartV0Glauber = 18.4951;
  Float_t errNpartCl2Glauber = 7.57; // 2% not used
  Float_t errNpartCl2Hij = 7.306;

  Float_t dNdEtaCorr = ncorrtrackletscentr/nevtscentr;
  Float_t errdNdEtaCorr = hCorrecteddNdEtaCentr->GetBinError(1);

  Float_t dNdEtaPartPairsV0G = dNdEtaCorr/(.5*NpartV0Glauber);
  Float_t dNdEtaPartPairsCl2G = dNdEtaCorr/(.5*NpartCl2Glauber);
  Float_t dNdEtaPartPairsCl2H = dNdEtaCorr/(.5*NpartCl2Hij);
  
  Float_t errdNdEtaPartPairsV0G  = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartV0Glauber,errdNdEtaCorr,errNpartV0Glauber);
  Float_t errdNdEtaPartPairsCl2G = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartCl2Glauber,errdNdEtaCorr,errNpartCl2Glauber);
  Float_t errdNdEtaPartPairsCl2H = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartCl2Hij,errdNdEtaCorr,errNpartCl2Hij);

  cout<<"V0 Glauber dN/dEta/.5<Npart> "<<dNdEtaPartPairsV0G<<" +- "<<errdNdEtaPartPairsV0G<<endl;
  cout<<"Cl2 Glauber dN/dEta/.5<Npart> "<<dNdEtaPartPairsCl2G<<" +- "<<errdNdEtaPartPairsCl2G<<endl;
  cout<<"Cl2 HIJING dN/dEta/.5<Npart> "<<dNdEtaPartPairsCl2H<<" +- "<<errdNdEtaPartPairsCl2H<<endl;

  hdNdEta_bkgCorr->Scale(nBinsPerPseudorapidityUnit);
  hdNdEta_fullyCorr->Scale(nBinsPerPseudorapidityUnit);

  // Create mask for MC prim in SPD acceptance
  TH1D *hdNdEta_test = new TH1D("dNdEta_test","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);  
  TH2F *hMask = new TH2F("Mask","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
    for (Int_t kvtx=0; kvtx<nBinVtx; kvtx++) {
       if (hSPDEta_2D_fullyCorr->GetBinContent(jeta+1,kvtx+1)!=0) hMask->SetBinContent(jeta+1,kvtx+1,1.);
    }
  }
  hMask->Multiply(Eta);
  hdNdEta_test->Add(hMask->ProjectionX("Mask_x",binVtxStart+1,binVtxStop,"e"));
  hdNdEta_test->Divide(hnorm_fullyCorr);
//  hdNdEta_test->Scale(1./nEvtsTotSel); // -> number of MC events  
  hdNdEta_test->Scale(nBinsPerPseudorapidityUnit);
  for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta_test->SetBinError(i+1,0.);

  TH1F* hRatiotest=new TH1F("ratiotest","",nBinEtaCorr,binLimEtaCorr);
  hRatiotest->GetXaxis()->SetTitle("#eta");
  hRatiotest->GetYaxis()->SetTitle("Generated/Corrected");
  hRatiotest->Divide(hdNdEta_test,hdNdEta_fullyCorr,1.,1.);
  new TCanvas();
  hRatiotest->Draw("p");

  // Generated/corrected ratios for consistency checks
  TH2F* hRatiodNdEta2D=new TH2F("ratiodNdEta2D","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hRatiodNdEta2D->Divide(hSPDEta_2D_fullyCorr,Eta,1.,1.);
  new TCanvas();
  hRatiodNdEta2D->Draw();

  TH1F* hRatiodNdEta=new TH1F("ratiodNdEta","",nBinEtaCorr,binLimEtaCorr);   
  hRatiodNdEta->GetXaxis()->SetTitle("#eta");
  hRatiodNdEta->GetYaxis()->SetTitle("Generated/Corrected");
  hRatiodNdEta->Divide(hdNdEta,hdNdEta_fullyCorr,1.,1.); 
  new TCanvas();
  hRatiodNdEta->Draw("p"); 

  // Draw dN/dEta
  new TCanvas();
//  hdNdEta->SetLineWidth(3);
//  hdNdEta->Draw("histo");

  hdNdEta_raw->SetMinimum(0.);
  hdNdEta_raw->SetMaximum(3000);
  hdNdEta_raw->SetLineColor(kGreen);
  hdNdEta_raw->SetLineWidth(3);
  hdNdEta_raw->Draw("histo");

  hdNdEta_bkgCorr->SetMarkerStyle(21);
  hdNdEta_bkgCorr->SetMarkerColor(kBlue);

  hdNdEta_bkgCorr->Draw("same,p");

  hdNdEta_fullyCorr->SetMarkerStyle(20);
  hdNdEta_fullyCorr->SetMarkerColor(kRed);
  hdNdEta_fullyCorr->Draw("p,same");

  TLegend *leg1 = new TLegend(0.2,0.7,0.7,0.9,NULL,"brNDC");
  leg1->SetFillColor(0);
  leg1->SetTextFont(62);
  leg1->SetTextSize(0.0304569);
  TLegendEntry *entry1=leg1->AddEntry(hdNdEta_raw,"Raw","l");
//                entry1=leg1->AddEntry(hdNdEta,"Generated","l");
                entry1=leg1->AddEntry(hdNdEta_bkgCorr,"Comb bkg corrected","p");
                entry1=leg1->AddEntry(hdNdEta_fullyCorr,"Eff/acc corrected","p");

  leg1->Draw();

/*  new TCanvas();
  // plot the relative stat error for this correction
  TH2F* hStatErrTrackToPart = new TH2F("staterrperctracktopart","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
  hStatErrTrackToPart->GetXaxis()->SetTitle("#eta");
  hStatErrTrackToPart->GetYaxis()->SetTitle("z_{vtx} [cm]");
  hStatErrTrackToPart->GetZaxis()->SetTitle("Statistical error (%)");
  for (Int_t kvtx=0; kvtx<nBinVtx; kvtx++) {
    for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
      if (hTrackToParticleCorr->GetBinContent(jeta+1,kvtx+1)) 
      hStatErrTrackToPart->SetBinContent(jeta+1,kvtx+1,(hTrackToParticleCorr->GetBinError(jeta+1,kvtx+1))/(hTrackToParticleCorr->GetBinContent(jeta+1,kvtx+1))*100);
      else hStatErrTrackToPart->SetBinContent(jeta+1,kvtx+1,0);
    }
  }
  new TCanvas();
  hStatErrTrackToPart->Draw();
  new TCanvas();
  hRatiodNdEta_trackToParticle->Draw();
  new TCanvas();
  hTrackToParticleCorr->Draw();*/

  // Write histos
  TFile *fout= new TFile("SPDdNdEta.root","RECREATE"); 

  hCombBkgCorrData->Write();
  hCombBkgCorrMC->Write();
  hCombBkgCorr2MC->Write();
  hAlphaCorr->Write();

  hdNdEta_raw->Write();
  hdNdEta_bkgCorr->Write();
  hdNdEta_fullyCorr->Write();  
 
  hCorrecteddNdEtaCentr->Write();

  hEffProtons->Write();
  hEffKaons->Write();
  hEffPions->Write();
  hEffOthers->Write();
 
  fout->Write();
  fout->Close();

}

Float_t CalculateErrordNdEtaPerPartPair(Float_t a, Float_t b, Float_t da, Float_t db) {

  Float_t errdndetaperpartpairs = (a/b)*TMath::Sqrt(TMath::Power(da,2)/TMath::Power(a,2)+TMath::Power(db,2)/TMath::Power(b,2)); 
  return errdndetaperpartpairs;

}

void combscalingfactors (TList *lall, TList *lcomb, TList *lmc, Double_t *scaleNormRange_rot, Double_t *scaleNormRange_labels) {


 TH2F *hall = (TH2F*) lall->FindObject("fHistSPDdePhideTheta");
 TH2F *hcomb = (TH2F*) lcomb->FindObject("fHistSPDdePhideTheta");

 TH2F *hlabel = (TH2F*) lmc->FindObject("fHistdPhidThetaComb");
 
 TH1D* hproall = new TH1D("_xall","",2000,-1.,1.);
 hproall->Sumw2();
 hproall = hall->ProjectionX("_xa",0,-1,"e");

 TH1D* hprocomb = new TH1D("_xcomb","",2000,-1.,1.);
 hprocomb->Sumw2();
 hprocomb = hcomb->ProjectionX("_xc",0,-1,"e");

 TH1D *hprolabel = new TH1D("_xcomblabel","",2000,-1.,1.);
 hprolabel = hlabel->ProjectionX("_xcombl",0,-1,"e");

 // Normalize to integral in [0.08,0.20]
 Double_t denNormRange = hproall->Integral(801,920) + hproall->Integral(1081,1200);
 Double_t numNormRange_rot   = hprocomb->Integral(801,920) + hprocomb->Integral(1081,1200);
 *scaleNormRange_rot = denNormRange / numNormRange_rot;

 new TCanvas();
 hproall->Draw("histo");

 hprocomb->Scale(*scaleNormRange_rot);
 cout<<"Scaling factor normalizing to close tails: "<<*scaleNormRange_rot<<endl;
 cout<<"Percentage of background"<< hprocomb->Integral(1,2000)/hproall->Integral(1,2000)<<endl;
 cout<<"Percentage of background in |#D#phi|<0.08 "<< hprocomb->Integral(921,1080)/hproall->Integral(921,1080)<<endl;
 cout<<"Percentage of background in |#D#phi|<0.06 "<< hprocomb->Integral(941,1060)/hproall->Integral(941,1060)<<endl;
 hprocomb->SetLineColor(kBlue);
 hprocomb->Draw("histo,same");

 // Fit the label bkg comb distribution to the data distribution
 // Normalize to integral in [0.08,0.20]
 Double_t numNormRange_labels   = hprolabel->Integral(801,920) + hprolabel->Integral(1081,1200);
 *scaleNormRange_labels = denNormRange / numNormRange_labels;

 hprolabel->Scale(*scaleNormRange_labels);
 cout<<"Scaling factor labels: "<<*scaleNormRange_labels<<endl;
 cout<<"Percentage of background with labels "<< (hprolabel->Integral(1,2000)/hproall->Integral(1,2000))<<endl;
 cout<<"Percentage of background with labels in |#D#phi|<0.08 "<< (hprolabel->Integral(921,1080)/hproall->Integral(921,1080))<<endl;

 hprolabel->SetLineColor(kGreen);
 hproall->Draw("histo,same");
/*
 TFile *fout= new TFile("scaling.root","RECREATE");
 hproall->Write();
 hprocomb->Write();
 hprolabel->Write();
 fout->Write();
 fout->Close();
*/
}


 correctSPDdNdEta.C:1
 correctSPDdNdEta.C:2
 correctSPDdNdEta.C:3
 correctSPDdNdEta.C:4
 correctSPDdNdEta.C:5
 correctSPDdNdEta.C:6
 correctSPDdNdEta.C:7
 correctSPDdNdEta.C:8
 correctSPDdNdEta.C:9
 correctSPDdNdEta.C:10
 correctSPDdNdEta.C:11
 correctSPDdNdEta.C:12
 correctSPDdNdEta.C:13
 correctSPDdNdEta.C:14
 correctSPDdNdEta.C:15
 correctSPDdNdEta.C:16
 correctSPDdNdEta.C:17
 correctSPDdNdEta.C:18
 correctSPDdNdEta.C:19
 correctSPDdNdEta.C:20
 correctSPDdNdEta.C:21
 correctSPDdNdEta.C:22
 correctSPDdNdEta.C:23
 correctSPDdNdEta.C:24
 correctSPDdNdEta.C:25
 correctSPDdNdEta.C:26
 correctSPDdNdEta.C:27
 correctSPDdNdEta.C:28
 correctSPDdNdEta.C:29
 correctSPDdNdEta.C:30
 correctSPDdNdEta.C:31
 correctSPDdNdEta.C:32
 correctSPDdNdEta.C:33
 correctSPDdNdEta.C:34
 correctSPDdNdEta.C:35
 correctSPDdNdEta.C:36
 correctSPDdNdEta.C:37
 correctSPDdNdEta.C:38
 correctSPDdNdEta.C:39
 correctSPDdNdEta.C:40
 correctSPDdNdEta.C:41
 correctSPDdNdEta.C:42
 correctSPDdNdEta.C:43
 correctSPDdNdEta.C:44
 correctSPDdNdEta.C:45
 correctSPDdNdEta.C:46
 correctSPDdNdEta.C:47
 correctSPDdNdEta.C:48
 correctSPDdNdEta.C:49
 correctSPDdNdEta.C:50
 correctSPDdNdEta.C:51
 correctSPDdNdEta.C:52
 correctSPDdNdEta.C:53
 correctSPDdNdEta.C:54
 correctSPDdNdEta.C:55
 correctSPDdNdEta.C:56
 correctSPDdNdEta.C:57
 correctSPDdNdEta.C:58
 correctSPDdNdEta.C:59
 correctSPDdNdEta.C:60
 correctSPDdNdEta.C:61
 correctSPDdNdEta.C:62
 correctSPDdNdEta.C:63
 correctSPDdNdEta.C:64
 correctSPDdNdEta.C:65
 correctSPDdNdEta.C:66
 correctSPDdNdEta.C:67
 correctSPDdNdEta.C:68
 correctSPDdNdEta.C:69
 correctSPDdNdEta.C:70
 correctSPDdNdEta.C:71
 correctSPDdNdEta.C:72
 correctSPDdNdEta.C:73
 correctSPDdNdEta.C:74
 correctSPDdNdEta.C:75
 correctSPDdNdEta.C:76
 correctSPDdNdEta.C:77
 correctSPDdNdEta.C:78
 correctSPDdNdEta.C:79
 correctSPDdNdEta.C:80
 correctSPDdNdEta.C:81
 correctSPDdNdEta.C:82
 correctSPDdNdEta.C:83
 correctSPDdNdEta.C:84
 correctSPDdNdEta.C:85
 correctSPDdNdEta.C:86
 correctSPDdNdEta.C:87
 correctSPDdNdEta.C:88
 correctSPDdNdEta.C:89
 correctSPDdNdEta.C:90
 correctSPDdNdEta.C:91
 correctSPDdNdEta.C:92
 correctSPDdNdEta.C:93
 correctSPDdNdEta.C:94
 correctSPDdNdEta.C:95
 correctSPDdNdEta.C:96
 correctSPDdNdEta.C:97
 correctSPDdNdEta.C:98
 correctSPDdNdEta.C:99
 correctSPDdNdEta.C:100
 correctSPDdNdEta.C:101
 correctSPDdNdEta.C:102
 correctSPDdNdEta.C:103
 correctSPDdNdEta.C:104
 correctSPDdNdEta.C:105
 correctSPDdNdEta.C:106
 correctSPDdNdEta.C:107
 correctSPDdNdEta.C:108
 correctSPDdNdEta.C:109
 correctSPDdNdEta.C:110
 correctSPDdNdEta.C:111
 correctSPDdNdEta.C:112
 correctSPDdNdEta.C:113
 correctSPDdNdEta.C:114
 correctSPDdNdEta.C:115
 correctSPDdNdEta.C:116
 correctSPDdNdEta.C:117
 correctSPDdNdEta.C:118
 correctSPDdNdEta.C:119
 correctSPDdNdEta.C:120
 correctSPDdNdEta.C:121
 correctSPDdNdEta.C:122
 correctSPDdNdEta.C:123
 correctSPDdNdEta.C:124
 correctSPDdNdEta.C:125
 correctSPDdNdEta.C:126
 correctSPDdNdEta.C:127
 correctSPDdNdEta.C:128
 correctSPDdNdEta.C:129
 correctSPDdNdEta.C:130
 correctSPDdNdEta.C:131
 correctSPDdNdEta.C:132
 correctSPDdNdEta.C:133
 correctSPDdNdEta.C:134
 correctSPDdNdEta.C:135
 correctSPDdNdEta.C:136
 correctSPDdNdEta.C:137
 correctSPDdNdEta.C:138
 correctSPDdNdEta.C:139
 correctSPDdNdEta.C:140
 correctSPDdNdEta.C:141
 correctSPDdNdEta.C:142
 correctSPDdNdEta.C:143
 correctSPDdNdEta.C:144
 correctSPDdNdEta.C:145
 correctSPDdNdEta.C:146
 correctSPDdNdEta.C:147
 correctSPDdNdEta.C:148
 correctSPDdNdEta.C:149
 correctSPDdNdEta.C:150
 correctSPDdNdEta.C:151
 correctSPDdNdEta.C:152
 correctSPDdNdEta.C:153
 correctSPDdNdEta.C:154
 correctSPDdNdEta.C:155
 correctSPDdNdEta.C:156
 correctSPDdNdEta.C:157
 correctSPDdNdEta.C:158
 correctSPDdNdEta.C:159
 correctSPDdNdEta.C:160
 correctSPDdNdEta.C:161
 correctSPDdNdEta.C:162
 correctSPDdNdEta.C:163
 correctSPDdNdEta.C:164
 correctSPDdNdEta.C:165
 correctSPDdNdEta.C:166
 correctSPDdNdEta.C:167
 correctSPDdNdEta.C:168
 correctSPDdNdEta.C:169
 correctSPDdNdEta.C:170
 correctSPDdNdEta.C:171
 correctSPDdNdEta.C:172
 correctSPDdNdEta.C:173
 correctSPDdNdEta.C:174
 correctSPDdNdEta.C:175
 correctSPDdNdEta.C:176
 correctSPDdNdEta.C:177
 correctSPDdNdEta.C:178
 correctSPDdNdEta.C:179
 correctSPDdNdEta.C:180
 correctSPDdNdEta.C:181
 correctSPDdNdEta.C:182
 correctSPDdNdEta.C:183
 correctSPDdNdEta.C:184
 correctSPDdNdEta.C:185
 correctSPDdNdEta.C:186
 correctSPDdNdEta.C:187
 correctSPDdNdEta.C:188
 correctSPDdNdEta.C:189
 correctSPDdNdEta.C:190
 correctSPDdNdEta.C:191
 correctSPDdNdEta.C:192
 correctSPDdNdEta.C:193
 correctSPDdNdEta.C:194
 correctSPDdNdEta.C:195
 correctSPDdNdEta.C:196
 correctSPDdNdEta.C:197
 correctSPDdNdEta.C:198
 correctSPDdNdEta.C:199
 correctSPDdNdEta.C:200
 correctSPDdNdEta.C:201
 correctSPDdNdEta.C:202
 correctSPDdNdEta.C:203
 correctSPDdNdEta.C:204
 correctSPDdNdEta.C:205
 correctSPDdNdEta.C:206
 correctSPDdNdEta.C:207
 correctSPDdNdEta.C:208
 correctSPDdNdEta.C:209
 correctSPDdNdEta.C:210
 correctSPDdNdEta.C:211
 correctSPDdNdEta.C:212
 correctSPDdNdEta.C:213
 correctSPDdNdEta.C:214
 correctSPDdNdEta.C:215
 correctSPDdNdEta.C:216
 correctSPDdNdEta.C:217
 correctSPDdNdEta.C:218
 correctSPDdNdEta.C:219
 correctSPDdNdEta.C:220
 correctSPDdNdEta.C:221
 correctSPDdNdEta.C:222
 correctSPDdNdEta.C:223
 correctSPDdNdEta.C:224
 correctSPDdNdEta.C:225
 correctSPDdNdEta.C:226
 correctSPDdNdEta.C:227
 correctSPDdNdEta.C:228
 correctSPDdNdEta.C:229
 correctSPDdNdEta.C:230
 correctSPDdNdEta.C:231
 correctSPDdNdEta.C:232
 correctSPDdNdEta.C:233
 correctSPDdNdEta.C:234
 correctSPDdNdEta.C:235
 correctSPDdNdEta.C:236
 correctSPDdNdEta.C:237
 correctSPDdNdEta.C:238
 correctSPDdNdEta.C:239
 correctSPDdNdEta.C:240
 correctSPDdNdEta.C:241
 correctSPDdNdEta.C:242
 correctSPDdNdEta.C:243
 correctSPDdNdEta.C:244
 correctSPDdNdEta.C:245
 correctSPDdNdEta.C:246
 correctSPDdNdEta.C:247
 correctSPDdNdEta.C:248
 correctSPDdNdEta.C:249
 correctSPDdNdEta.C:250
 correctSPDdNdEta.C:251
 correctSPDdNdEta.C:252
 correctSPDdNdEta.C:253
 correctSPDdNdEta.C:254
 correctSPDdNdEta.C:255
 correctSPDdNdEta.C:256
 correctSPDdNdEta.C:257
 correctSPDdNdEta.C:258
 correctSPDdNdEta.C:259
 correctSPDdNdEta.C:260
 correctSPDdNdEta.C:261
 correctSPDdNdEta.C:262
 correctSPDdNdEta.C:263
 correctSPDdNdEta.C:264
 correctSPDdNdEta.C:265
 correctSPDdNdEta.C:266
 correctSPDdNdEta.C:267
 correctSPDdNdEta.C:268
 correctSPDdNdEta.C:269
 correctSPDdNdEta.C:270
 correctSPDdNdEta.C:271
 correctSPDdNdEta.C:272
 correctSPDdNdEta.C:273
 correctSPDdNdEta.C:274
 correctSPDdNdEta.C:275
 correctSPDdNdEta.C:276
 correctSPDdNdEta.C:277
 correctSPDdNdEta.C:278
 correctSPDdNdEta.C:279
 correctSPDdNdEta.C:280
 correctSPDdNdEta.C:281
 correctSPDdNdEta.C:282
 correctSPDdNdEta.C:283
 correctSPDdNdEta.C:284
 correctSPDdNdEta.C:285
 correctSPDdNdEta.C:286
 correctSPDdNdEta.C:287
 correctSPDdNdEta.C:288
 correctSPDdNdEta.C:289
 correctSPDdNdEta.C:290
 correctSPDdNdEta.C:291
 correctSPDdNdEta.C:292
 correctSPDdNdEta.C:293
 correctSPDdNdEta.C:294
 correctSPDdNdEta.C:295
 correctSPDdNdEta.C:296
 correctSPDdNdEta.C:297
 correctSPDdNdEta.C:298
 correctSPDdNdEta.C:299
 correctSPDdNdEta.C:300
 correctSPDdNdEta.C:301
 correctSPDdNdEta.C:302
 correctSPDdNdEta.C:303
 correctSPDdNdEta.C:304
 correctSPDdNdEta.C:305
 correctSPDdNdEta.C:306
 correctSPDdNdEta.C:307
 correctSPDdNdEta.C:308
 correctSPDdNdEta.C:309
 correctSPDdNdEta.C:310
 correctSPDdNdEta.C:311
 correctSPDdNdEta.C:312
 correctSPDdNdEta.C:313
 correctSPDdNdEta.C:314
 correctSPDdNdEta.C:315
 correctSPDdNdEta.C:316
 correctSPDdNdEta.C:317
 correctSPDdNdEta.C:318
 correctSPDdNdEta.C:319
 correctSPDdNdEta.C:320
 correctSPDdNdEta.C:321
 correctSPDdNdEta.C:322
 correctSPDdNdEta.C:323
 correctSPDdNdEta.C:324
 correctSPDdNdEta.C:325
 correctSPDdNdEta.C:326
 correctSPDdNdEta.C:327
 correctSPDdNdEta.C:328
 correctSPDdNdEta.C:329
 correctSPDdNdEta.C:330
 correctSPDdNdEta.C:331
 correctSPDdNdEta.C:332
 correctSPDdNdEta.C:333
 correctSPDdNdEta.C:334
 correctSPDdNdEta.C:335
 correctSPDdNdEta.C:336
 correctSPDdNdEta.C:337
 correctSPDdNdEta.C:338
 correctSPDdNdEta.C:339
 correctSPDdNdEta.C:340
 correctSPDdNdEta.C:341
 correctSPDdNdEta.C:342
 correctSPDdNdEta.C:343
 correctSPDdNdEta.C:344
 correctSPDdNdEta.C:345
 correctSPDdNdEta.C:346
 correctSPDdNdEta.C:347
 correctSPDdNdEta.C:348
 correctSPDdNdEta.C:349
 correctSPDdNdEta.C:350
 correctSPDdNdEta.C:351
 correctSPDdNdEta.C:352
 correctSPDdNdEta.C:353
 correctSPDdNdEta.C:354
 correctSPDdNdEta.C:355
 correctSPDdNdEta.C:356
 correctSPDdNdEta.C:357
 correctSPDdNdEta.C:358
 correctSPDdNdEta.C:359
 correctSPDdNdEta.C:360
 correctSPDdNdEta.C:361
 correctSPDdNdEta.C:362
 correctSPDdNdEta.C:363
 correctSPDdNdEta.C:364
 correctSPDdNdEta.C:365
 correctSPDdNdEta.C:366
 correctSPDdNdEta.C:367
 correctSPDdNdEta.C:368
 correctSPDdNdEta.C:369
 correctSPDdNdEta.C:370
 correctSPDdNdEta.C:371
 correctSPDdNdEta.C:372
 correctSPDdNdEta.C:373
 correctSPDdNdEta.C:374
 correctSPDdNdEta.C:375
 correctSPDdNdEta.C:376
 correctSPDdNdEta.C:377
 correctSPDdNdEta.C:378
 correctSPDdNdEta.C:379
 correctSPDdNdEta.C:380
 correctSPDdNdEta.C:381
 correctSPDdNdEta.C:382
 correctSPDdNdEta.C:383
 correctSPDdNdEta.C:384
 correctSPDdNdEta.C:385
 correctSPDdNdEta.C:386
 correctSPDdNdEta.C:387
 correctSPDdNdEta.C:388
 correctSPDdNdEta.C:389
 correctSPDdNdEta.C:390
 correctSPDdNdEta.C:391
 correctSPDdNdEta.C:392
 correctSPDdNdEta.C:393
 correctSPDdNdEta.C:394
 correctSPDdNdEta.C:395
 correctSPDdNdEta.C:396
 correctSPDdNdEta.C:397
 correctSPDdNdEta.C:398
 correctSPDdNdEta.C:399
 correctSPDdNdEta.C:400
 correctSPDdNdEta.C:401
 correctSPDdNdEta.C:402
 correctSPDdNdEta.C:403
 correctSPDdNdEta.C:404
 correctSPDdNdEta.C:405
 correctSPDdNdEta.C:406
 correctSPDdNdEta.C:407
 correctSPDdNdEta.C:408
 correctSPDdNdEta.C:409
 correctSPDdNdEta.C:410
 correctSPDdNdEta.C:411
 correctSPDdNdEta.C:412
 correctSPDdNdEta.C:413
 correctSPDdNdEta.C:414
 correctSPDdNdEta.C:415
 correctSPDdNdEta.C:416
 correctSPDdNdEta.C:417
 correctSPDdNdEta.C:418
 correctSPDdNdEta.C:419
 correctSPDdNdEta.C:420
 correctSPDdNdEta.C:421
 correctSPDdNdEta.C:422
 correctSPDdNdEta.C:423
 correctSPDdNdEta.C:424
 correctSPDdNdEta.C:425
 correctSPDdNdEta.C:426
 correctSPDdNdEta.C:427
 correctSPDdNdEta.C:428
 correctSPDdNdEta.C:429
 correctSPDdNdEta.C:430
 correctSPDdNdEta.C:431
 correctSPDdNdEta.C:432
 correctSPDdNdEta.C:433
 correctSPDdNdEta.C:434
 correctSPDdNdEta.C:435
 correctSPDdNdEta.C:436
 correctSPDdNdEta.C:437
 correctSPDdNdEta.C:438
 correctSPDdNdEta.C:439
 correctSPDdNdEta.C:440
 correctSPDdNdEta.C:441
 correctSPDdNdEta.C:442
 correctSPDdNdEta.C:443
 correctSPDdNdEta.C:444
 correctSPDdNdEta.C:445
 correctSPDdNdEta.C:446
 correctSPDdNdEta.C:447
 correctSPDdNdEta.C:448
 correctSPDdNdEta.C:449
 correctSPDdNdEta.C:450
 correctSPDdNdEta.C:451
 correctSPDdNdEta.C:452
 correctSPDdNdEta.C:453
 correctSPDdNdEta.C:454
 correctSPDdNdEta.C:455
 correctSPDdNdEta.C:456
 correctSPDdNdEta.C:457
 correctSPDdNdEta.C:458
 correctSPDdNdEta.C:459
 correctSPDdNdEta.C:460
 correctSPDdNdEta.C:461
 correctSPDdNdEta.C:462
 correctSPDdNdEta.C:463
 correctSPDdNdEta.C:464
 correctSPDdNdEta.C:465
 correctSPDdNdEta.C:466
 correctSPDdNdEta.C:467
 correctSPDdNdEta.C:468
 correctSPDdNdEta.C:469
 correctSPDdNdEta.C:470
 correctSPDdNdEta.C:471
 correctSPDdNdEta.C:472
 correctSPDdNdEta.C:473
 correctSPDdNdEta.C:474
 correctSPDdNdEta.C:475
 correctSPDdNdEta.C:476
 correctSPDdNdEta.C:477
 correctSPDdNdEta.C:478
 correctSPDdNdEta.C:479
 correctSPDdNdEta.C:480
 correctSPDdNdEta.C:481
 correctSPDdNdEta.C:482
 correctSPDdNdEta.C:483
 correctSPDdNdEta.C:484
 correctSPDdNdEta.C:485
 correctSPDdNdEta.C:486
 correctSPDdNdEta.C:487
 correctSPDdNdEta.C:488
 correctSPDdNdEta.C:489
 correctSPDdNdEta.C:490
 correctSPDdNdEta.C:491
 correctSPDdNdEta.C:492
 correctSPDdNdEta.C:493
 correctSPDdNdEta.C:494
 correctSPDdNdEta.C:495
 correctSPDdNdEta.C:496
 correctSPDdNdEta.C:497
 correctSPDdNdEta.C:498
 correctSPDdNdEta.C:499
 correctSPDdNdEta.C:500
 correctSPDdNdEta.C:501
 correctSPDdNdEta.C:502
 correctSPDdNdEta.C:503
 correctSPDdNdEta.C:504
 correctSPDdNdEta.C:505
 correctSPDdNdEta.C:506
 correctSPDdNdEta.C:507
 correctSPDdNdEta.C:508
 correctSPDdNdEta.C:509
 correctSPDdNdEta.C:510
 correctSPDdNdEta.C:511
 correctSPDdNdEta.C:512
 correctSPDdNdEta.C:513
 correctSPDdNdEta.C:514
 correctSPDdNdEta.C:515
 correctSPDdNdEta.C:516
 correctSPDdNdEta.C:517
 correctSPDdNdEta.C:518
 correctSPDdNdEta.C:519
 correctSPDdNdEta.C:520
 correctSPDdNdEta.C:521
 correctSPDdNdEta.C:522
 correctSPDdNdEta.C:523
 correctSPDdNdEta.C:524
 correctSPDdNdEta.C:525
 correctSPDdNdEta.C:526
 correctSPDdNdEta.C:527
 correctSPDdNdEta.C:528
 correctSPDdNdEta.C:529
 correctSPDdNdEta.C:530
 correctSPDdNdEta.C:531
 correctSPDdNdEta.C:532
 correctSPDdNdEta.C:533
 correctSPDdNdEta.C:534
 correctSPDdNdEta.C:535
 correctSPDdNdEta.C:536
 correctSPDdNdEta.C:537
 correctSPDdNdEta.C:538
 correctSPDdNdEta.C:539
 correctSPDdNdEta.C:540
 correctSPDdNdEta.C:541
 correctSPDdNdEta.C:542
 correctSPDdNdEta.C:543
 correctSPDdNdEta.C:544
 correctSPDdNdEta.C:545
 correctSPDdNdEta.C:546
 correctSPDdNdEta.C:547
 correctSPDdNdEta.C:548
 correctSPDdNdEta.C:549
 correctSPDdNdEta.C:550
 correctSPDdNdEta.C:551
 correctSPDdNdEta.C:552
 correctSPDdNdEta.C:553
 correctSPDdNdEta.C:554
 correctSPDdNdEta.C:555
 correctSPDdNdEta.C:556
 correctSPDdNdEta.C:557
 correctSPDdNdEta.C:558
 correctSPDdNdEta.C:559
 correctSPDdNdEta.C:560
 correctSPDdNdEta.C:561
 correctSPDdNdEta.C:562
 correctSPDdNdEta.C:563
 correctSPDdNdEta.C:564
 correctSPDdNdEta.C:565
 correctSPDdNdEta.C:566
 correctSPDdNdEta.C:567
 correctSPDdNdEta.C:568
 correctSPDdNdEta.C:569
 correctSPDdNdEta.C:570
 correctSPDdNdEta.C:571
 correctSPDdNdEta.C:572
 correctSPDdNdEta.C:573
 correctSPDdNdEta.C:574
 correctSPDdNdEta.C:575
 correctSPDdNdEta.C:576
 correctSPDdNdEta.C:577
 correctSPDdNdEta.C:578
 correctSPDdNdEta.C:579
 correctSPDdNdEta.C:580
 correctSPDdNdEta.C:581
 correctSPDdNdEta.C:582
 correctSPDdNdEta.C:583
 correctSPDdNdEta.C:584
 correctSPDdNdEta.C:585
 correctSPDdNdEta.C:586
 correctSPDdNdEta.C:587
 correctSPDdNdEta.C:588
 correctSPDdNdEta.C:589
 correctSPDdNdEta.C:590
 correctSPDdNdEta.C:591
 correctSPDdNdEta.C:592
 correctSPDdNdEta.C:593
 correctSPDdNdEta.C:594
 correctSPDdNdEta.C:595