[u/mrichter/AliRoot.git] / PWGLF / totEt / macros / calculateCorrections.C

//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 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)
{
  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 *gammaEff2D;
  //if(fHistGammasGeneratedMult && fHistGammasFoundMult){
    gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,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 );
  

  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;
	numeratorVal /= rescale;
      }
      Double_t denominatorValErr = denominator->GetBinError(ibin,jbin);
      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,jbin,quotient);
      result->SetBinError(ibin,jbin,quotientErr);
      //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
    }
  }
  return result;
}
Commit	Line	Data
9f3689ba	1	//Oystein Djuvsland
9f3689ba	2	//macro for plotting correction factors for EM et for use with AliAnalysisEt and writes them to a file
943dc129	3	#include "TTree.h"
	4	#include "TFile.h"
	5	#include <TList.h>
	6	#include <Rtypes.h>
	7	#include "TCanvas.h"
	8	#include <TH2I.h>
	9	#include <TH2F.h>
	10	#include <TProfile.h>
	11	#include <TFitResultPtr.h>
	12	#include <TFitResult.h>
	13	#include "TF1.h"
	14	#include <iostream>
	15	#include <AliAnalysisEtTrackMatchCorrections.h>
	16	#include <AliAnalysisEtRecEffCorrection.h>
	17
	18	TCanvas *c1 = 0;
	19
47151f26	20	TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ;
47151f26	21	TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ;
9f3689ba	22
9f3689ba	23	//creates an empty set of correction factors for debugging purposes
943dc129	24	TF1 * generateRecEffFunction(Double_t p0, Double_t p1);
9f3689ba	25	int createDummy(Double_t p0 = 1.0, Double_t p1 = 0.0, char *det = "Phos")
02c62614	26	{
	27	TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
	28	TF1 fitneutral("fitneutral","0", 0, 100);
	29	TF1 fitcharged("fitcharged","0", 0, 100);
	30	TF1 fitgamma("fitgamma","0", 0, 100);
	31	TF1 fitsecondary("fitsecondary","0", 0, 100);
9f3689ba	32	AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(Form("TmCorrections%s",det), fitcharged, fitneutral, fitgamma, fitsecondary,0,0,0,0);
02c62614	33	TF1 *func = generateRecEffFunction(p0, p1);
9f3689ba	34	AliAnalysisEtRecEffCorrection recor = new AliAnalysisEtRecEffCorrection(Form("ReCorrections%s",det), func, 1000);
02c62614	35
	36	cor->Write();
	37	recor->Write();
	38	outfile->Close();
	39
943dc129	40
02c62614	41	}
9f3689ba	42	//p0 = correction factors for efficiency from track matching
	43	//p1 = correction factors for efficiency from track matching
	44	//determined from fit of efficiency from track matching, 0.366 from simulations, fit as a function of energy
	45	//p0 is the constant
	46	//p1 is linear term
9793e78d	47	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)
943dc129	48	{
9f3689ba	49	TString detector = det;
9793e78d	50	TString emcal = "Emcal";
943dc129	51	c1 = new TCanvas;
	52
	53	TFile *f = TFile::Open(filename, "READ");
	54
	55	TList l = (TList)f->Get("out1");
	56
9f3689ba	57	TTree primaryTree = (TTree)l->FindObject(("fPrimaryTree"+detector+"MC").Data());
	58	TTree recTree = (TTree)l->FindObject(("fEventSummaryTree"+detector+"Rec").Data());
	59	TTree mcTree = (TTree)l->FindObject(("fEventSummaryTree"+detector+"MC").Data());
943dc129	60	std::cout << primaryTree << " " << recTree << " " << mcTree << std::endl;
	61
	62	Int_t clusterMult = 0;
	63	Int_t nChargedNotRemoved = 0;
	64	Int_t nNeutralNotRemoved = 0;
	65	Int_t nGammaRemoved = 0;
	66	Int_t nSecNotRemoved = 0;
	67
	68	Double_t emEtRec = 0;
	69	Double_t emEtMc = 0;
	70
	71	recTree->SetBranchAddress("fNeutralMultiplicity", &clusterMult);
	72	mcTree->SetBranchAddress("fChargedNotRemoved", &nChargedNotRemoved);
	73	mcTree->SetBranchAddress("fNeutralNotRemoved", &nNeutralNotRemoved);
	74	mcTree->SetBranchAddress("fGammaRemoved", &nGammaRemoved);
	75	mcTree->SetBranchAddress("fSecondaryNotRemoved", &nSecNotRemoved);
	76
9793e78d	77	Float_t maxMult = 99.5;
9793e78d	78	Int_t nbins = 100;
d45a33b8	79	if(detector==emcal){
9793e78d	80	// //100 seems to be sufficient...
d45a33b8	81	maxMult = 249.5;
	82	nbins = 250;
	83	}
9793e78d	84	TH2I *hChargedVsClusterMult = new TH2I("hChVsMult", "hChVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
	85	TH2I *hNeutralVsClusterMult = new TH2I("hNeutVsMult", "hNeutVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
	86	TH2I *hGammaVsClusterMult = new TH2I("hGammaVsMult", "hGammaVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
	87	TH2I *hSecVsClusterMult = new TH2I("hSecVsMult", "hSecVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
	88
943dc129	89	Int_t nEvents = mcTree->GetEntriesFast();
	90	for(Int_t i = 0; i < nEvents; i++)
	91	{
	92	mcTree->GetEvent(i);
	93	recTree->GetEvent(i);
	94	hChargedVsClusterMult->Fill(clusterMult, nChargedNotRemoved);
	95	hNeutralVsClusterMult->Fill(clusterMult, nNeutralNotRemoved);
	96	hGammaVsClusterMult->Fill(clusterMult, nGammaRemoved);
	97	hSecVsClusterMult->Fill(clusterMult, nSecNotRemoved);
	98	}
	99
	100	c1->Divide(2,2);
	101	c1->cd(1);
	102	TString title = "Charged particles not removed";
	103	TString xtitle = "Cluster Multiplicity";
	104	TString ytitle = "N_{ch}";
	105	hChargedVsClusterMult->SetTitle(title);
	106	hChargedVsClusterMult->GetYaxis()->SetTitle(ytitle);
	107	hChargedVsClusterMult->GetXaxis()->SetTitle(xtitle);
	108	hChargedVsClusterMult->SetStats(0);
	109	hChargedVsClusterMult->Draw();
	110	TProfile *chProf = hChargedVsClusterMult->ProfileX();
	111	chProf->SetStats(0);
	112	chProf->Draw("SAME");
	113	//TF1 fitcharged("fitcharged","([0] + [1]x)(0.48/([2] + [3]*[2]))", 0, 100);
9793e78d	114	TF1 fitcharged("fitcharged","pol2", 0, 100);//fit of number of charged tracks vs detector multiplicity
9f3689ba	115	//if straight line track matching roughly not dependent on centrality
943dc129	116	// fitcharged.FixParameter(2, p0);
	117	// fitcharged.FixParameter(3, p1);
	118	TFitResultPtr chRes = chProf->Fit(&fitcharged,"S");
	119	TArrayD ch;
	120	if(!chRes)
	121	{
	122	ch = TArrayD(chRes->NPar(),chRes->GetParams());
	123	}
	124	else
	125	{
	126	std::cout << "Could not extract charged contribution params" << std::endl;
	127	}
	128	c1->cd(2);
	129	title = "Neutral particles not removed";
	130	ytitle = "N_{neutral}";
	131	hNeutralVsClusterMult->SetTitle(title);
	132	hNeutralVsClusterMult->GetXaxis()->SetTitle(xtitle);
	133	hNeutralVsClusterMult->GetYaxis()->SetTitle(ytitle);
	134	hNeutralVsClusterMult->SetStats(0);
	135	hNeutralVsClusterMult->Draw();
	136	TProfile *neuProf = hNeutralVsClusterMult->ProfileX();
	137	neuProf->SetStats(0);
	138	neuProf->Draw("SAME");
9793e78d	139	TF1 fitneutral("fitneutral","pol2", 0, 100);//fit of number of neutral particles in calo that we call background in calo vs detector multiplicity
9f3689ba	140	//may include K0S
943dc129	141	TFitResultPtr neuRes = neuProf->Fit(&fitneutral,"S");
	142	TArrayD neu;
	143	if(!neuRes)
	144	{
	145	neu = TArrayD(neuRes->NPar(),neuRes->GetParams());
	146	}
	147	else
	148	{
	149	std::cout << "Could not extract charged contribution params" << std::endl;
	150	}
	151	c1->cd(3);
	152
	153	title = "Gammas removed";
	154	ytitle = "N_{#gamma}";
	155	hGammaVsClusterMult->SetTitle(title);
	156	hGammaVsClusterMult->GetYaxis()->SetTitle(ytitle);
	157	hGammaVsClusterMult->GetXaxis()->SetTitle(xtitle);
	158	hGammaVsClusterMult->SetStats(0);
	159	hGammaVsClusterMult->Draw();
	160	TProfile *gamProf = hGammaVsClusterMult->ProfileX();
	161	gamProf->SetStats(0);
	162	gamProf->Draw("SAME");
9793e78d	163	TF1 fitgamma("fitgamma","pol2", 0, 100);//fit of number of gammas removed erroneously vs detector multiplicity
943dc129	164	TFitResultPtr gammaRes = gamProf->Fit(&fitgamma,"S");
	165	TArrayD gamma;
	166	if(!gammaRes)
	167	{
	168	gamma = TArrayD(gammaRes->NPar(),gammaRes->GetParams());
	169	}
	170	else
	171	{
	172	std::cout << "Could not extract charged contribution params" << std::endl;
	173	}
	174
	175	c1->cd(4);
	176
	177	title = "Secondaries not removed";
	178	ytitle = "N_{sec}";
	179	hSecVsClusterMult->SetTitle(title);
	180	hSecVsClusterMult->GetXaxis()->SetTitle(xtitle);
	181	hSecVsClusterMult->GetYaxis()->SetTitle(ytitle);
	182	hSecVsClusterMult->SetStats(0);
	183	hSecVsClusterMult->Draw();
	184	TProfile *secProf = hSecVsClusterMult->ProfileX();
	185	secProf->SetStats(0);
	186	secProf->Draw("SAME");
9793e78d	187	TF1 fitsecondary("fitsecondary","pol2", 0, 100);//fit of number of secondary particles that leave deposits in calo vs detector multiplicity
943dc129	188	TFitResultPtr secRes = secProf->Fit(&fitsecondary,"S");
	189	TArrayD sec;
	190	if(!secRes)
	191	{
	192	sec = TArrayD(secRes->NPar(),secRes->GetParams());
	193	}
	194	else
	195	{
	196	std::cout << "Could not extract charged contribution params" << std::endl;
	197	}
9f3689ba	198	//ugly hack for getting the energy
	199	//changing number of particles to energy of particles by multiplying by the <energy> and dividing by the efficiency for the average energy
	200	//average energy from each of these: in excel file
943dc129	201	Double_t meanCharged = 0.48/(p0 + p1*0.48);
	202	Double_t meanNeutral = 0.53/(p0 + p1*0.53);
	203	Double_t meanGamma = 0.51/(p0 + p1*0.51);
	204	Double_t meanSecondary = meanGamma;
	205
95309578	206	TH2F *fHistHadronDepositsAllMult = l->FindObject("fHistHadronDepositsAllCent");
	207	TH2F *fHistHadronDepositsRecoMult = l->FindObject("fHistHadronDepositsRecoCent");
	208	TH2F *eff2D;
d45a33b8	209	// if(fHistHadronDepositsRecoMult && fHistHadronDepositsAllMult ){
95309578	210	eff2D = (TH2F*) bayneseffdiv2D(fHistHadronDepositsRecoMult,fHistHadronDepositsAllMult,"eff2D");
d45a33b8	211	// }
	212	// else{
	213	// cerr<<"Warning! Did not calculate reconstruction efficiency!!"<<endl;
	214	// eff2D = new TH2F("eff2D", "eff2D",200, 0, 10,20,-0.5,19.5);
	215	// }
47151f26	216	//cor->SetReconstructionEfficiency(eff2D);
47151f26	217
95309578	218	TH2F *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedCent");
	219	TH2F *fHistGammasFoundMult = l->FindObject("fHistGammasFoundCent");
	220	TH2F *gammaEff2D;
d45a33b8	221	//if(fHistGammasGeneratedMult && fHistGammasFoundMult){
95309578	222	gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,fHistGammasGeneratedMult,"gammaEff2D");
d45a33b8	223	// }
	224	// else{
	225	// cerr<<"Warning! Did not calculate reconstruction efficiency!!"<<endl;
	226	// gammaEff2D = new TH2F("gammaEff2D", "gammaEff2D",200, 0, 10,20,-0.5,19.5);
	227	// }
47151f26	228
47151f26	229	AliAnalysisEtTrackMatchCorrections cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, eff2D,
943dc129	230	meanCharged, meanNeutral, meanGamma, meanSecondary );
943dc129	231
47151f26	232
943dc129	233	TF1 *func = generateRecEffFunction(p0, p1);
47151f26	234	AliAnalysisEtRecEffCorrection recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), func,*gammaEff2D, 1000);
943dc129	235	TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
	236	cor->Write();
	237	recor->Write();
	238	return 0;
	239	}
	240
	241
	242	TF1* generateRecEffFunction(Double_t p0, Double_t p1)
	243	{
	244	TF1 f = new TF1("receff", "x/([0] + x[1])", 0, 200);
	245
	246	Double_t params[2] = {p0, p1};
	247	f->SetParameters(params);
	248	return f;
	249	}
47151f26	250
	251	TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
	252	{
	253	if(!numerator){
	254	cerr<<"Error: numerator does not exist!"<<endl;
	255	return NULL;
	256	}
	257	if(!denominator){
	258	cerr<<"Error: denominator does not exist!"<<endl;
	259	return NULL;
	260	}
	261	TH1* result = (TH1*)numerator->Clone(name);
	262	Int_t nbins = numerator->GetNbinsX();
	263	for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
	264	Double_t numeratorVal = numerator->GetBinContent(ibin);
	265	Double_t denominatorVal = denominator->GetBinContent(ibin);
	266	// Check if the errors are right or the thing is scaled
	267	Double_t numeratorValErr = numerator->GetBinError(ibin);
	268	if (!(numeratorValErr==0. \|\| numeratorVal ==0.) ) {
	269	Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
	270	numeratorVal /= rescale;
	271	}
	272	Double_t denominatorValErr = denominator->GetBinError(ibin);
	273	if (!(denominatorValErr==0. \|\| denominatorVal==0. )) {
	274	Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
	275	denominatorVal /= rescale;
	276	}
	277	Double_t quotient = 0.;
	278	if (denominatorVal!=0.) {
	279	quotient = numeratorVal/denominatorVal;
	280	}
	281	Double_t quotientErr=0;
	282	quotientErr = TMath::Sqrt(
	283	(numeratorVal+1.0)/(denominatorVal+2.0)*
	284	((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
	285	result->SetBinContent(ibin,quotient);
	286	result->SetBinError(ibin,quotientErr);
	287	//cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
	288	}
	289	return result;
	290	}
	291
	292
	293	TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name)
	294	{
	295	if(!numerator){
	296	cerr<<"Error: numerator does not exist!"<<endl;
	297	return NULL;
	298	}
	299	if(!denominator){
	300	cerr<<"Error: denominator does not exist!"<<endl;
	301	return NULL;
	302	}
	303	TH2* result = (TH2*)numerator->Clone(name);
	304	Int_t nbinsX = numerator->GetNbinsX();
	305	Int_t nbinsY = numerator->GetNbinsY();
	306	for (Int_t ibin=0; ibin<= nbinsX+1; ++ibin) {
	307	for (Int_t jbin=0; jbin<= nbinsY+1; ++jbin) {
	308	Double_t numeratorVal = numerator->GetBinContent(ibin,jbin);
	309	Double_t denominatorVal = denominator->GetBinContent(ibin,jbin);
	310	// Check if the errors are right or the thing is scaled
	311	Double_t numeratorValErr = numerator->GetBinError(ibin,jbin);
	312	if (!(numeratorValErr==0. \|\| numeratorVal ==0.) ) {
	313	Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
314	numeratorVal /= rescale;
315	}
316	Double_t denominatorValErr = denominator->GetBinError(ibin,jbin);
317	if (!(denominatorValErr==0. \|\| denominatorVal==0. )) {
318	Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
319	denominatorVal /= rescale;
320	}
321	Double_t quotient = 0.;
322	if (denominatorVal!=0.) {
323	quotient = numeratorVal/denominatorVal;
324	}
325	Double_t quotientErr=0;
326	quotientErr = TMath::Sqrt(
327	(numeratorVal+1.0)/(denominatorVal+2.0)*
328	((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
329	result->SetBinContent(ibin,jbin,quotient);
330	result->SetBinError(ibin,jbin,quotientErr);
331	//cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
332	}
333	}
334	return result;
335	}