ROOT logo
//Oystein Djuvsland
//macro for plotting correction factors for EM et for use with AliAnalysisEt and writes them to a file
#include "TTree.h"
#include "TFile.h"
#include <TList.h>
#include <Rtypes.h>
#include "TCanvas.h"
#include <TH2I.h>
#include <TH2F.h>
#include <TProfile.h>
#include <TFitResultPtr.h>
#include <TFitResult.h>
#include "TF1.h"
#include <iostream>
#include <AliAnalysisEtTrackMatchCorrections.h>
#include <AliAnalysisEtRecEffCorrection.h>

TCanvas *c1 = 0;

TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ;
TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ;

//creates an empty set of correction factors for debugging purposes
TF1 * generateRecEffFunction(Double_t p0, Double_t p1);
int createDummy(Double_t p0 = 1.0, Double_t p1 = 0.0, char *det = "Phos")
{
  TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
  TF1 fitneutral("fitneutral","0", 0, 100);
  TF1 fitcharged("fitcharged","0", 0, 100);
  TF1 fitgamma("fitgamma","0", 0, 100);
  TF1 fitsecondary("fitsecondary","0", 0, 100);
  AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(Form("TmCorrections%s",det), fitcharged, fitneutral, fitgamma, fitsecondary,0,0,0,0);
  TF1 *func = generateRecEffFunction(p0, p1);
  AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(Form("ReCorrections%s",det), *func, 1000);

  cor->Write();
  recor->Write();
  outfile->Close();


}
//p0 = correction factors for efficiency from track matching
//p1 = correction factors for efficiency from track matching
//determined from fit of efficiency from track matching, 0.366 from simulations, fit as a function of energy
//p0 is the constant
//p1 is linear term
Int_t CutSet = 0;//Defaults: 250 MeV for PHOS and 300 MeV for EMCal
int calculateCorrections(TString filename="rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root", Double_t p0 = 0.366, Double_t p1 = 0.0, char *det = "Emcal", Bool_t isSim = kFALSE)
{
  string inline;
  float value = 0;
  float error = 0;
  int i=0;
  TString detector = det;
  TString emcal = "Emcal";
  c1 = new TCanvas;
  
  TFile *f = TFile::Open(filename, "READ");
  
  TList *l = (TList*)f->Get("out1");
  
//   TTree *primaryTree = (TTree*)l->FindObject(("fPrimaryTree"+detector+"MC").Data());
//   TTree *recTree = (TTree*)l->FindObject(("fEventSummaryTree"+detector+"Rec").Data());
//   TTree *mcTree = (TTree*)l->FindObject(("fEventSummaryTree"+detector+"MC").Data());
//   std::cout << primaryTree << " " << recTree << " " << mcTree << std::endl;
 
  Int_t clusterMult = 0;
  Int_t nChargedNotRemoved = 0;
  Int_t nNeutralNotRemoved = 0;
  Int_t nGammaRemoved = 0;
  Int_t nSecNotRemoved = 0;
  
  Double_t emEtRec = 0;
  Double_t emEtMc = 0;

//   recTree->SetBranchAddress("fNeutralMultiplicity", &clusterMult);
//   mcTree->SetBranchAddress("fChargedNotRemoved", &nChargedNotRemoved);
//   mcTree->SetBranchAddress("fNeutralNotRemoved", &nNeutralNotRemoved);
//   mcTree->SetBranchAddress("fGammaRemoved", &nGammaRemoved);
//   mcTree->SetBranchAddress("fSecondaryNotRemoved", &nSecNotRemoved);
  
  Float_t maxMult = 99.5;
  Int_t nbins = 100;
   if(detector==emcal){
//     //100 seems to be sufficient...
     maxMult = 249.5;
     nbins = 250;
  }
  TH2I *hChargedVsClusterMult = new TH2I("hChVsMult", "hChVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
  TH2I *hNeutralVsClusterMult = new TH2I("hNeutVsMult", "hNeutVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
  TH2I *hGammaVsClusterMult = new TH2I("hGammaVsMult", "hGammaVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
  TH2I *hSecVsClusterMult = new TH2I("hSecVsMult", "hSecVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);

//   Int_t nEvents = mcTree->GetEntriesFast();
//     for(Int_t i = 0; i < nEvents; i++)
//     {
//       mcTree->GetEvent(i);
//       recTree->GetEvent(i);
//       hChargedVsClusterMult->Fill(clusterMult, nChargedNotRemoved);
//       hNeutralVsClusterMult->Fill(clusterMult, nNeutralNotRemoved);
//       hGammaVsClusterMult->Fill(clusterMult, nGammaRemoved);
//       hSecVsClusterMult->Fill(clusterMult, nSecNotRemoved);
//     }
    
  c1->Divide(2,2);
  c1->cd(1);
  TString title = "Charged particles not removed";
  TString xtitle = "Cluster Multiplicity";
  TString ytitle = "N_{ch}";
  hChargedVsClusterMult->SetTitle(title);
  hChargedVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hChargedVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hChargedVsClusterMult->SetStats(0);
  hChargedVsClusterMult->Draw();
  TProfile *chProf = hChargedVsClusterMult->ProfileX();
  chProf->SetStats(0);
  chProf->Draw("SAME");
  //TF1 fitcharged("fitcharged","([0] + [1]*x)*(0.48/([2] + [3]*[2]))", 0, 100);
  TF1 fitcharged("fitcharged","pol2", 0, 100);//fit of number of charged tracks vs detector multiplicity
  //if straight line track matching roughly not dependent on centrality
  //fitcharged.FixParameter(2, p0);
   //   fitcharged.FixParameter(3, p1);
  TFitResultPtr chRes = chProf->Fit(&fitcharged,"S");
  TArrayD ch;
  if(!chRes)
  {
    ch = TArrayD(chRes->NPar(),chRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  c1->cd(2);
  title = "Neutral particles not removed";
  ytitle = "N_{neutral}";
  hNeutralVsClusterMult->SetTitle(title);
  hNeutralVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hNeutralVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hNeutralVsClusterMult->SetStats(0);
  hNeutralVsClusterMult->Draw();
  TProfile *neuProf = hNeutralVsClusterMult->ProfileX();
  neuProf->SetStats(0);
  neuProf->Draw("SAME");
  TF1 fitneutral("fitneutral","pol2", 0, 100);//fit of number of neutral particles in calo that we call background in calo vs detector multiplicity
  //may include K0S
  TFitResultPtr neuRes = neuProf->Fit(&fitneutral,"S");
  TArrayD neu;
  if(!neuRes)
  {
    neu = TArrayD(neuRes->NPar(),neuRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  c1->cd(3);
  
  title = "Gammas removed";
  ytitle = "N_{#gamma}";
  hGammaVsClusterMult->SetTitle(title);
  hGammaVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hGammaVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hGammaVsClusterMult->SetStats(0);
  hGammaVsClusterMult->Draw();
  TProfile *gamProf = hGammaVsClusterMult->ProfileX();
  gamProf->SetStats(0);
  gamProf->Draw("SAME");
  TF1 fitgamma("fitgamma","pol2", 0, 100);//fit of number of gammas removed erroneously vs detector multiplicity
  TFitResultPtr gammaRes = gamProf->Fit(&fitgamma,"S");
  TArrayD gamma;
  if(!gammaRes)
  {
    gamma = TArrayD(gammaRes->NPar(),gammaRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  
  c1->cd(4);
  
  title = "Secondaries not removed";
  ytitle = "N_{sec}";
  hSecVsClusterMult->SetTitle(title);
  hSecVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hSecVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hSecVsClusterMult->SetStats(0);
  hSecVsClusterMult->Draw();
  TProfile *secProf = hSecVsClusterMult->ProfileX();
  secProf->SetStats(0);
  secProf->Draw("SAME");
  TF1 fitsecondary("fitsecondary","pol2", 0, 100);//fit of number of secondary particles that leave deposits in calo vs detector multiplicity
  TFitResultPtr secRes = secProf->Fit(&fitsecondary,"S");
  TArrayD sec;
  if(!secRes)
  {
    sec = TArrayD(secRes->NPar(),secRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  //ugly hack for getting the energy
  //changing number of particles to energy of particles by multiplying by the <energy> and dividing by the efficiency for the average energy
  //average energy from each of these:  in excel file
  Double_t meanCharged = 0.48/(p0 + p1*0.48);
  Double_t meanNeutral = 0.53/(p0 + p1*0.53);
  Double_t meanGamma = 0.51/(p0 + p1*0.51);
  Double_t meanSecondary = meanGamma; 
  
  TH2F  *fHistHadronDepositsAllMult = l->FindObject("fHistHadronDepositsAllCent");
  TH2F  *fHistHadronDepositsRecoMult = l->FindObject("fHistHadronDepositsRecoCent");
  TH2F *eff2D;
//   if(fHistHadronDepositsRecoMult && fHistHadronDepositsAllMult ){
    eff2D = (TH2F*) bayneseffdiv2D(fHistHadronDepositsRecoMult,fHistHadronDepositsAllMult,"eff2D");
//   }
//   else{
//     cerr<<"Warning!  Did not calculate reconstruction efficiency!!"<<endl;
//     eff2D =  new TH2F("eff2D", "eff2D",200, 0, 10,20,-0.5,19.5);
//   }
  //cor->SetReconstructionEfficiency(eff2D);
  
  TH2F  *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedCent");
//   TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundCent");
//   TH2F  *fHistGammasFoundOutOfAccMult = l->FindObject("fHistGammasFoundOutOfAccCent");
   TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundRecoEnergyCent");
   TH2F  *fHistGammasFoundOutOfAccMult = l->FindObject("fHistGammasFoundOutOfAccRecoEnergyCent");
  if(fHistGammasFoundOutOfAccMult){
    //cout<<"I have "<<fHistGammasFoundOutOfAccMult->GetEntries()<<" entries"<<endl;
    fHistGammasFoundMult->Add(fHistGammasFoundOutOfAccMult);
    //cout<<"I have "<<fHistGammasFoundMult->GetEntries()<<" entries"<<endl;
  }
  else{cout<<"fHistGammasFoundOutOfAccCent does not exist!"<<endl;}
  //fHistGammasFoundMult->Add(fHistGammasFoundOutOfAccMult);
  TH2F *gammaEff2D;
  //if(fHistGammasGeneratedMult && fHistGammasFoundMult){
  gammaEff2D = (TH2F*) bayneseffdiv2D((TH2F*)fHistGammasFoundMult,(TH2F*)fHistGammasGeneratedMult,"gammaEff2D");
//   }
//   else{
//     cerr<<"Warning!  Did not calculate reconstruction efficiency!!"<<endl;
//     gammaEff2D =  new TH2F("gammaEff2D", "gammaEff2D",200, 0, 10,20,-0.5,19.5);
//   }

  AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, *eff2D,
										   meanCharged, meanNeutral, meanGamma, meanSecondary );
  
  TString corrfilename = "allcorrections"+detector+".dat";
  cout<<"Reading "<<corrfilename<<endl;
  ifstream myfile (corrfilename.Data());
  Float_t neutroncorr = 0;
  Float_t hadroncorr = 0;
  Float_t kaoncorr = 0;
  Float_t secondarycorr = 0;
  Float_t minetcorr = 0;
  Float_t junk = 0;
  i=0;
  if (myfile.is_open()){
    while ( myfile.good() )
      {
	getline (myfile,inline);
	istringstream tmp(inline);
// 	tmp >> neutroncorr;
// 	tmp >> hadroncorr;
// 	tmp >> kaoncorr;
// 	tmp >> secondarycorr;

	tmp >> junk;
	tmp >> junk;
	tmp >> junk;
	tmp >> junk;
	tmp >> minetcorr;
	tmp >> junk;//phosMinEtError[i];
	tmp >> junk;//phosNonLinError[i];
	tmp >> neutroncorr;//phosNeutronCorr[i];
	tmp >> junk;//phosNeutronError[i];
	tmp >> hadroncorr;//phosHadronCorr[i];
	tmp >> junk;//phosHadronError[i];
	tmp >> kaoncorr;//phosKaonCorr[i];
	tmp >> junk;//phosKaonError[i];
	tmp >> secondarycorr;//phosSecondaryCorr[i];
	tmp >> junk;//phosSecondaryError[i];
	if(i<20){
	  //cout<<" cb "<<i<<" "<<minetcorr<<" "<<neutroncorr<<" "<<hadroncorr<<" "<<kaoncorr<<" "<<secondarycorr<<endl;
	  cor->SetMinEtCorrection(i,minetcorr);
 	  cor->SetNeutronCorrection(i,neutroncorr);
 	  cor->SetHadronCorrection(i,hadroncorr);
 	  cor->SetKaonCorrection(i,kaoncorr);
 	  cor->SetSecondaryCorrection(i,secondarycorr);
	  i++;

	}
      }
    myfile.close();
  }


//   TString cutstring = "";
//   if(CutSet==1) cutstring = "Cut6";
//   if(CutSet==2) cutstring = "Cut7";
//   TString minetInfileNameShort = "MinEt"+detector+"Short"+cutstring+".dat";
//   cout<<"Reading "<<minetInfileNameShort.Data()<<endl;
//   Float_t minEtErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
//   Float_t minEtCorrShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
//   ifstream myminetfileShort (minetInfileNameShort.Data());
//   if (myminetfileShort.is_open()){
//     while ( myminetfileShort.good() )
//       {
// 	getline (myminetfileShort,inline);
// 	istringstream tmp(inline);
// 	tmp >> value;
// 	tmp >> error;
// 	if(i<10){
// 	  minEtCorrShort[i] = value;
// 	  minEtErrorShort[i] = error;
// 	}
// 	cout<<"min et corr cb "<<i<<" "<< value <<" +/- "<< error <<endl;
// 	i++;
//       }
//     myminetfileShort.close();
//   }



  TF1 *func = generateRecEffFunction(p0, p1);
  AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), *func,*gammaEff2D, 1000);
  TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
  cor->Write();
  recor->Write();
  return 0;
}


TF1* generateRecEffFunction(Double_t p0, Double_t p1)
{
  TF1 *f = new TF1("receff", "x/([0] + x*[1])", 0, 200);
  
  Double_t params[2] = {p0, p1};
  f->SetParameters(params);
  return f;
}

TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) 
{
    if(!numerator){
      cerr<<"Error:  numerator does not exist!"<<endl;
      return NULL;
    }
    if(!denominator){
      cerr<<"Error:  denominator does not exist!"<<endl;
      return NULL;
    }
    TH1* result = (TH1*)numerator->Clone(name);
    Int_t nbins = numerator->GetNbinsX();
    for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
      Double_t numeratorVal = numerator->GetBinContent(ibin);
      Double_t denominatorVal = denominator->GetBinContent(ibin);
      // Check if the errors are right or the thing is scaled
      Double_t numeratorValErr = numerator->GetBinError(ibin);
      if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
	Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
	numeratorVal /= rescale;
      }
      Double_t denominatorValErr = denominator->GetBinError(ibin);
      if (!(denominatorValErr==0. || denominatorVal==0. )) {
	Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
	denominatorVal /= rescale;
      }
      Double_t quotient = 0.;
      if (denominatorVal!=0.) {
	quotient = numeratorVal/denominatorVal;
      }
      Double_t quotientErr=0;
      quotientErr = TMath::Sqrt(
				(numeratorVal+1.0)/(denominatorVal+2.0)*
				((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
      result->SetBinContent(ibin,quotient);
      result->SetBinError(ibin,quotientErr);
      //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
    }
    return result;
}


TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) 
{
  if(!numerator){
    cerr<<"Error:  numerator does not exist!"<<endl;
    return NULL;
  }
  if(!denominator){
    cerr<<"Error:  denominator does not exist!"<<endl;
    return NULL;
  }
  TH2* result = (TH2*)numerator->Clone(name);
  Int_t nbinsX = numerator->GetNbinsX();
  Int_t nbinsY = numerator->GetNbinsY();
  for (Int_t ibin=0; ibin<= nbinsX+1; ++ibin) {
    for (Int_t jbin=0; jbin<= nbinsY+1; ++jbin) {
      Double_t numeratorVal = numerator->GetBinContent(ibin,jbin);
      Double_t denominatorVal = denominator->GetBinContent(ibin,jbin);
      // Check if the errors are right or the thing is scaled
      Double_t numeratorValErr = numerator->GetBinError(ibin,jbin);
      if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
	Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
	if(rescale != 0.0) numeratorVal /= rescale;
      }
      Double_t denominatorValErr = denominator->GetBinError(ibin,jbin);
      if (!(denominatorValErr==0. || denominatorVal==0. )) {
	Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
	if(rescale != 0.0) denominatorVal /= rescale;
      }
      Double_t quotient = 0.;
      if (denominatorVal!=0.) {
	quotient = numeratorVal/denominatorVal;
      }
      Double_t quotientErr=0;
      if(denominatorVal>0){
	//cerr<<(numeratorVal+1.0)<<"/"<<(denominatorVal+2.0)<<"*"<<"("<<(numeratorVal+2.0)<<"/"<<(denominatorVal+3.0)<<"-"<<(numeratorVal+1.0)<<"/"<<(denominatorVal+2.0)<<")"<<endl;
	float val = (numeratorVal+1.0)/(denominatorVal+2.0)*
	  ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0));
	if(val>0) quotientErr = TMath::Sqrt(val);
      }
      result->SetBinContent(ibin,jbin,quotient);
      result->SetBinError(ibin,jbin,quotientErr);
      //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
    }
  }
  return result;
}
 calculateCorrections.C:1
 calculateCorrections.C:2
 calculateCorrections.C:3
 calculateCorrections.C:4
 calculateCorrections.C:5
 calculateCorrections.C:6
 calculateCorrections.C:7
 calculateCorrections.C:8
 calculateCorrections.C:9
 calculateCorrections.C:10
 calculateCorrections.C:11
 calculateCorrections.C:12
 calculateCorrections.C:13
 calculateCorrections.C:14
 calculateCorrections.C:15
 calculateCorrections.C:16
 calculateCorrections.C:17
 calculateCorrections.C:18
 calculateCorrections.C:19
 calculateCorrections.C:20
 calculateCorrections.C:21
 calculateCorrections.C:22
 calculateCorrections.C:23
 calculateCorrections.C:24
 calculateCorrections.C:25
 calculateCorrections.C:26
 calculateCorrections.C:27
 calculateCorrections.C:28
 calculateCorrections.C:29
 calculateCorrections.C:30
 calculateCorrections.C:31
 calculateCorrections.C:32
 calculateCorrections.C:33
 calculateCorrections.C:34
 calculateCorrections.C:35
 calculateCorrections.C:36
 calculateCorrections.C:37
 calculateCorrections.C:38
 calculateCorrections.C:39
 calculateCorrections.C:40
 calculateCorrections.C:41
 calculateCorrections.C:42
 calculateCorrections.C:43
 calculateCorrections.C:44
 calculateCorrections.C:45
 calculateCorrections.C:46
 calculateCorrections.C:47
 calculateCorrections.C:48
 calculateCorrections.C:49
 calculateCorrections.C:50
 calculateCorrections.C:51
 calculateCorrections.C:52
 calculateCorrections.C:53
 calculateCorrections.C:54
 calculateCorrections.C:55
 calculateCorrections.C:56
 calculateCorrections.C:57
 calculateCorrections.C:58
 calculateCorrections.C:59
 calculateCorrections.C:60
 calculateCorrections.C:61
 calculateCorrections.C:62
 calculateCorrections.C:63
 calculateCorrections.C:64
 calculateCorrections.C:65
 calculateCorrections.C:66
 calculateCorrections.C:67
 calculateCorrections.C:68
 calculateCorrections.C:69
 calculateCorrections.C:70
 calculateCorrections.C:71
 calculateCorrections.C:72
 calculateCorrections.C:73
 calculateCorrections.C:74
 calculateCorrections.C:75
 calculateCorrections.C:76
 calculateCorrections.C:77
 calculateCorrections.C:78
 calculateCorrections.C:79
 calculateCorrections.C:80
 calculateCorrections.C:81
 calculateCorrections.C:82
 calculateCorrections.C:83
 calculateCorrections.C:84
 calculateCorrections.C:85
 calculateCorrections.C:86
 calculateCorrections.C:87
 calculateCorrections.C:88
 calculateCorrections.C:89
 calculateCorrections.C:90
 calculateCorrections.C:91
 calculateCorrections.C:92
 calculateCorrections.C:93
 calculateCorrections.C:94
 calculateCorrections.C:95
 calculateCorrections.C:96
 calculateCorrections.C:97
 calculateCorrections.C:98
 calculateCorrections.C:99
 calculateCorrections.C:100
 calculateCorrections.C:101
 calculateCorrections.C:102
 calculateCorrections.C:103
 calculateCorrections.C:104
 calculateCorrections.C:105
 calculateCorrections.C:106
 calculateCorrections.C:107
 calculateCorrections.C:108
 calculateCorrections.C:109
 calculateCorrections.C:110
 calculateCorrections.C:111
 calculateCorrections.C:112
 calculateCorrections.C:113
 calculateCorrections.C:114
 calculateCorrections.C:115
 calculateCorrections.C:116
 calculateCorrections.C:117
 calculateCorrections.C:118
 calculateCorrections.C:119
 calculateCorrections.C:120
 calculateCorrections.C:121
 calculateCorrections.C:122
 calculateCorrections.C:123
 calculateCorrections.C:124
 calculateCorrections.C:125
 calculateCorrections.C:126
 calculateCorrections.C:127
 calculateCorrections.C:128
 calculateCorrections.C:129
 calculateCorrections.C:130
 calculateCorrections.C:131
 calculateCorrections.C:132
 calculateCorrections.C:133
 calculateCorrections.C:134
 calculateCorrections.C:135
 calculateCorrections.C:136
 calculateCorrections.C:137
 calculateCorrections.C:138
 calculateCorrections.C:139
 calculateCorrections.C:140
 calculateCorrections.C:141
 calculateCorrections.C:142
 calculateCorrections.C:143
 calculateCorrections.C:144
 calculateCorrections.C:145
 calculateCorrections.C:146
 calculateCorrections.C:147
 calculateCorrections.C:148
 calculateCorrections.C:149
 calculateCorrections.C:150
 calculateCorrections.C:151
 calculateCorrections.C:152
 calculateCorrections.C:153
 calculateCorrections.C:154
 calculateCorrections.C:155
 calculateCorrections.C:156
 calculateCorrections.C:157
 calculateCorrections.C:158
 calculateCorrections.C:159
 calculateCorrections.C:160
 calculateCorrections.C:161
 calculateCorrections.C:162
 calculateCorrections.C:163
 calculateCorrections.C:164
 calculateCorrections.C:165
 calculateCorrections.C:166
 calculateCorrections.C:167
 calculateCorrections.C:168
 calculateCorrections.C:169
 calculateCorrections.C:170
 calculateCorrections.C:171
 calculateCorrections.C:172
 calculateCorrections.C:173
 calculateCorrections.C:174
 calculateCorrections.C:175
 calculateCorrections.C:176
 calculateCorrections.C:177
 calculateCorrections.C:178
 calculateCorrections.C:179
 calculateCorrections.C:180
 calculateCorrections.C:181
 calculateCorrections.C:182
 calculateCorrections.C:183
 calculateCorrections.C:184
 calculateCorrections.C:185
 calculateCorrections.C:186
 calculateCorrections.C:187
 calculateCorrections.C:188
 calculateCorrections.C:189
 calculateCorrections.C:190
 calculateCorrections.C:191
 calculateCorrections.C:192
 calculateCorrections.C:193
 calculateCorrections.C:194
 calculateCorrections.C:195
 calculateCorrections.C:196
 calculateCorrections.C:197
 calculateCorrections.C:198
 calculateCorrections.C:199
 calculateCorrections.C:200
 calculateCorrections.C:201
 calculateCorrections.C:202
 calculateCorrections.C:203
 calculateCorrections.C:204
 calculateCorrections.C:205
 calculateCorrections.C:206
 calculateCorrections.C:207
 calculateCorrections.C:208
 calculateCorrections.C:209
 calculateCorrections.C:210
 calculateCorrections.C:211
 calculateCorrections.C:212
 calculateCorrections.C:213
 calculateCorrections.C:214
 calculateCorrections.C:215
 calculateCorrections.C:216
 calculateCorrections.C:217
 calculateCorrections.C:218
 calculateCorrections.C:219
 calculateCorrections.C:220
 calculateCorrections.C:221
 calculateCorrections.C:222
 calculateCorrections.C:223
 calculateCorrections.C:224
 calculateCorrections.C:225
 calculateCorrections.C:226
 calculateCorrections.C:227
 calculateCorrections.C:228
 calculateCorrections.C:229
 calculateCorrections.C:230
 calculateCorrections.C:231
 calculateCorrections.C:232
 calculateCorrections.C:233
 calculateCorrections.C:234
 calculateCorrections.C:235
 calculateCorrections.C:236
 calculateCorrections.C:237
 calculateCorrections.C:238
 calculateCorrections.C:239
 calculateCorrections.C:240
 calculateCorrections.C:241
 calculateCorrections.C:242
 calculateCorrections.C:243
 calculateCorrections.C:244
 calculateCorrections.C:245
 calculateCorrections.C:246
 calculateCorrections.C:247
 calculateCorrections.C:248
 calculateCorrections.C:249
 calculateCorrections.C:250
 calculateCorrections.C:251
 calculateCorrections.C:252
 calculateCorrections.C:253
 calculateCorrections.C:254
 calculateCorrections.C:255
 calculateCorrections.C:256
 calculateCorrections.C:257
 calculateCorrections.C:258
 calculateCorrections.C:259
 calculateCorrections.C:260
 calculateCorrections.C:261
 calculateCorrections.C:262
 calculateCorrections.C:263
 calculateCorrections.C:264
 calculateCorrections.C:265
 calculateCorrections.C:266
 calculateCorrections.C:267
 calculateCorrections.C:268
 calculateCorrections.C:269
 calculateCorrections.C:270
 calculateCorrections.C:271
 calculateCorrections.C:272
 calculateCorrections.C:273
 calculateCorrections.C:274
 calculateCorrections.C:275
 calculateCorrections.C:276
 calculateCorrections.C:277
 calculateCorrections.C:278
 calculateCorrections.C:279
 calculateCorrections.C:280
 calculateCorrections.C:281
 calculateCorrections.C:282
 calculateCorrections.C:283
 calculateCorrections.C:284
 calculateCorrections.C:285
 calculateCorrections.C:286
 calculateCorrections.C:287
 calculateCorrections.C:288
 calculateCorrections.C:289
 calculateCorrections.C:290
 calculateCorrections.C:291
 calculateCorrections.C:292
 calculateCorrections.C:293
 calculateCorrections.C:294
 calculateCorrections.C:295
 calculateCorrections.C:296
 calculateCorrections.C:297
 calculateCorrections.C:298
 calculateCorrections.C:299
 calculateCorrections.C:300
 calculateCorrections.C:301
 calculateCorrections.C:302
 calculateCorrections.C:303
 calculateCorrections.C:304
 calculateCorrections.C:305
 calculateCorrections.C:306
 calculateCorrections.C:307
 calculateCorrections.C:308
 calculateCorrections.C:309
 calculateCorrections.C:310
 calculateCorrections.C:311
 calculateCorrections.C:312
 calculateCorrections.C:313
 calculateCorrections.C:314
 calculateCorrections.C:315
 calculateCorrections.C:316
 calculateCorrections.C:317
 calculateCorrections.C:318
 calculateCorrections.C:319
 calculateCorrections.C:320
 calculateCorrections.C:321
 calculateCorrections.C:322
 calculateCorrections.C:323
 calculateCorrections.C:324
 calculateCorrections.C:325
 calculateCorrections.C:326
 calculateCorrections.C:327
 calculateCorrections.C:328
 calculateCorrections.C:329
 calculateCorrections.C:330
 calculateCorrections.C:331
 calculateCorrections.C:332
 calculateCorrections.C:333
 calculateCorrections.C:334
 calculateCorrections.C:335
 calculateCorrections.C:336
 calculateCorrections.C:337
 calculateCorrections.C:338
 calculateCorrections.C:339
 calculateCorrections.C:340
 calculateCorrections.C:341
 calculateCorrections.C:342
 calculateCorrections.C:343
 calculateCorrections.C:344
 calculateCorrections.C:345
 calculateCorrections.C:346
 calculateCorrections.C:347
 calculateCorrections.C:348
 calculateCorrections.C:349
 calculateCorrections.C:350
 calculateCorrections.C:351
 calculateCorrections.C:352
 calculateCorrections.C:353
 calculateCorrections.C:354
 calculateCorrections.C:355
 calculateCorrections.C:356
 calculateCorrections.C:357
 calculateCorrections.C:358
 calculateCorrections.C:359
 calculateCorrections.C:360
 calculateCorrections.C:361
 calculateCorrections.C:362
 calculateCorrections.C:363
 calculateCorrections.C:364
 calculateCorrections.C:365
 calculateCorrections.C:366
 calculateCorrections.C:367
 calculateCorrections.C:368
 calculateCorrections.C:369
 calculateCorrections.C:370
 calculateCorrections.C:371
 calculateCorrections.C:372
 calculateCorrections.C:373
 calculateCorrections.C:374
 calculateCorrections.C:375
 calculateCorrections.C:376
 calculateCorrections.C:377
 calculateCorrections.C:378
 calculateCorrections.C:379
 calculateCorrections.C:380
 calculateCorrections.C:381
 calculateCorrections.C:382
 calculateCorrections.C:383
 calculateCorrections.C:384
 calculateCorrections.C:385
 calculateCorrections.C:386
 calculateCorrections.C:387
 calculateCorrections.C:388
 calculateCorrections.C:389
 calculateCorrections.C:390
 calculateCorrections.C:391
 calculateCorrections.C:392
 calculateCorrections.C:393
 calculateCorrections.C:394
 calculateCorrections.C:395
 calculateCorrections.C:396
 calculateCorrections.C:397
 calculateCorrections.C:398
 calculateCorrections.C:399
 calculateCorrections.C:400
 calculateCorrections.C:401
 calculateCorrections.C:402
 calculateCorrections.C:403
 calculateCorrections.C:404
 calculateCorrections.C:405
 calculateCorrections.C:406
 calculateCorrections.C:407
 calculateCorrections.C:408
 calculateCorrections.C:409
 calculateCorrections.C:410
 calculateCorrections.C:411
 calculateCorrections.C:412
 calculateCorrections.C:413
 calculateCorrections.C:414
 calculateCorrections.C:415
 calculateCorrections.C:416
 calculateCorrections.C:417
 calculateCorrections.C:418
 calculateCorrections.C:419
 calculateCorrections.C:420
 calculateCorrections.C:421
 calculateCorrections.C:422
 calculateCorrections.C:423
 calculateCorrections.C:424
 calculateCorrections.C:425
 calculateCorrections.C:426
 calculateCorrections.C:427
 calculateCorrections.C:428
 calculateCorrections.C:429
 calculateCorrections.C:430