[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisHadEt.cxx

//Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
//University of Tennessee at Knoxville
//
// This class is designed for the analysis of the hadronic component of 
// transverse energy.  It is used by AliAnalysisTaskHadEt.
// This gets information about the hadronic component of the transverse energy 
// from tracks reconstructed in an event
// it has daughters, AliAnalysisHadEtMonteCarlo and 
// AliAnalysisHadEtReconstructed which loop over either Monte Carlo data or 
// real data to get Et

#include "AliAnalysisHadEt.h"
#include "TMath.h"
#include "TList.h"
#include "TH1F.h"
#include "TH2F.h"
#include <iostream>
#include "AliAnalysisEtCuts.h"
#include "AliMCEvent.h"
#include "AliVEvent.h"
#include "AliStack.h"
#include "AliESDtrackCuts.h"
#include "TDatabasePDG.h"
#include "TParticle.h"
#include "Rtypes.h"

using namespace std;

ClassImp(AliAnalysisHadEt);


Int_t AliAnalysisHadEt::fgnumOfEtaBins = 46;
Float_t AliAnalysisHadEt::fgEtaAxis[47]={-0.78, -0.74, -0.7, -0.66, -0.62, -0.58, -0.54, -0.5, -0.46, -0.42, -0.38, -0.34, -0.3, -0.26, -0.22, -0.18, -0.14, -0.12, -0.1, -0.08, -0.06, -0.04, -0.02, -0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.18, 0.22, 0.26, 0.3, 0.34, 0.38, 0.42, 0.46, 0.5, 0.54, 0.58, 0.62, 0.66, 0.7, 0.74, 0.78};
Int_t AliAnalysisHadEt::fgNumOfPtBins = 111;
Float_t AliAnalysisHadEt::fgPtAxis[117]=
  {0.0,0.01,0.02,0.03,0.04, 0.05, 0.06,0.07,0.08,0.09, 0.10,0.11, .12,0.13, .14,0.15, .16,0.17, .18,0.19,
   0.2, .22, .24, .26, .28, 0.30, 0.32, .34, .36, .38, 0.40, .42, .44, .46, .48,
   0.5, .52, .54, .56, .58, 0.60, 0.62, .64, .66, .68, 0.70, .72, .74, .76, .78,
   .80, .82, .84, .86, .88, 0.90, 0.92, .94, .96, .98, 1.00,1.05, 1.1,1.15, 1.2,
  1.25, 1.3,1.35,1.40,1.45, 1.50, 1.55, 1.6,1.65, 1.7, 1.75, 1.8,1.85, 1.9,1.95,
   2.0, 2.2, 2.4, 2.6, 2.8, 3.00, 3.20, 3.4, 3.6, 3.8, 4.00, 4.2, 4.4, 4.6, 4.8,
   5.0, 5.5, 6.0, 6.5, 7.0, 7.50, 8.00, 8.5, 9.0, 9.5, 10.0,12.0,14.0,16.0,18.0,
  20.0,25.0,30.0,35.0,40.0, 45.0, 50.0}; 

Float_t AliAnalysisHadEt::fgPtTPCCutOff = 0.15;
Float_t AliAnalysisHadEt::fgPtITSCutOff = 0.10;

AliAnalysisHadEt::AliAnalysisHadEt() :
        fHistogramNameSuffix("")
	,fCuts(0)
        ,fPdgDB(0)
        ,fPiPlusCode(0)
        ,fPiMinusCode(0)
        ,fKPlusCode(0)
        ,fKMinusCode(0)
        ,fProtonCode(0)
        ,fAntiProtonCode(0)
        ,fLambdaCode(0)
        ,fAntiLambdaCode(0)
        ,fK0SCode(0)
        ,fOmegaCode(0)
        ,fAntiOmegaCode(0)
        ,fXi0Code(0)
        ,fAntiXi0Code(0)
        ,fXiCode(0)
        ,fAntiXiCode(0)
        ,fSigmaCode(0)
        ,fAntiSigmaCode(0)
        ,fK0LCode(0)
        ,fNeutronCode(0)
        ,fAntiNeutronCode(0)
        ,fEPlusCode(0)
        ,fEMinusCode(0)
	,fGammaCode(0)
	,fPi0Code(0)
	,fEtaCode(0)
	,fOmega0Code(0)
        ,fPionMass(0)
        ,fKaonMass(0)
        ,fProtonMass(0)
        ,fElectronMass(0)
	,fSumEt(0)
	,fSumEtAcc(0)
	,fTotEt(0)
	,fTotEtAcc(0)
	,fTotNeutralEt(0)
	,fTotNeutralEtAcc(0)
	,fTotChargedEt(0)
	,fTotChargedEtAcc(0)
	,fMultiplicity(0)
	,fChargedMultiplicity(0)
	,fNeutralMultiplicity(0)
        ,fEsdtrackCutsITSTPC(0)
        ,fEsdtrackCutsTPC(0)
	,fEsdtrackCutsITS(0)
        ,fhistoList(0)
{//default constructor

}

AliAnalysisHadEt::~AliAnalysisHadEt()
{//destructor

}

Int_t AliAnalysisHadEt::AnalyseEvent(AliVEvent *event)
{ //this line is basically here to eliminate a compiler warning that event is not used.  Making it a virtual function did not work with the plugin.
  cout << "This event has " << event->GetNumberOfTracks() << " tracks" << endl;
  ResetEventValues();
  return 0;
}

void AliAnalysisHadEt::FillOutputList()
{//fill the output histogram list with histograms in all AliAnalysisHadEt's
}

void AliAnalysisHadEt::Init()
{// clear variables, set up cuts and PDG info
  ResetEventValues();

}

void AliAnalysisHadEt::CreateHistograms()
{//creates histograms included in all AliAnalysisHadEt's
}

void AliAnalysisHadEt::FillHistograms()
{//Fills histograms filled for all AliAnalysisHadEt's
}

void AliAnalysisHadEt::ResetEventValues()
{//Resets event values of et to zero
  fTotEt = 0;
  fTotEtAcc = 0;
  fTotNeutralEt = 0;
  fTotNeutralEtAcc = 0;
  fTotChargedEt  = 0;
  fTotChargedEtAcc = 0;
  fMultiplicity = 0;
  fChargedMultiplicity = 0;
  fNeutralMultiplicity = 0;
  
  if (!fCuts || !fPdgDB || fPiPlusCode==0) { // some Init's needed
    cout << __FILE__ << ":" << __LINE__ << " : Init " << endl;
    if (!fCuts) {
      cout << " setting up Cuts " << endl;
      fCuts = new AliAnalysisEtCuts();
    }
    if(!fPdgDB) {
      cout << " setting up PdgDB " << endl;
      fPdgDB = new TDatabasePDG();
    }
    
    if (fPiPlusCode==0) {
      SetParticleCodes();
    }
  }
}

void AliAnalysisHadEt::SetParticleCodes()
{  //the codes are defined in $ROOTSYS/etc/pdg_table.txt
    fPionMass = fPdgDB->GetParticle("pi+")->Mass();
    fKaonMass = fPdgDB->GetParticle("K+")->Mass();
    fProtonMass = fPdgDB->GetParticle("proton")->Mass();
    fElectronMass = fPdgDB->GetParticle("e+")->Mass();
    fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
    fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
    fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
    fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
    fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
    fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
    fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
    fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
    fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
    fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
    fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
    fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
    fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
    fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
    fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
    fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
    fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
    fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
    fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
    fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
    fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
    fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
    fGammaCode =  fPdgDB->GetParticle("gamma")->PdgCode();
    fPi0Code =  fPdgDB->GetParticle("pi0")->PdgCode();
    fEtaCode =  fPdgDB->GetParticle("eta")->PdgCode();
    fOmega0Code =  fPdgDB->GetParticle("omega")->PdgCode();
    cout << "Resetting Codes: Pion " << fPiPlusCode
	 << "," << fPiMinusCode 
	 << " Kaon " << fKPlusCode 
	 << "," << fKMinusCode 
	 << " Gamma "<<fGammaCode
	 << " Pi0 "<<fPi0Code
	 << " Eta "<<fEtaCode
	 << " omega "<<fOmega0Code
	 << endl;
}

void AliAnalysisHadEt::CreateEtaPtHisto2D(TString name, TString title)
{     //creates a 2-d histogram in eta and phi and adds it to the list of histograms to be saved
  TString *histoname   = new TString();
  TString *histotitle   = new TString();

  histoname->Append(name);
  histotitle->Append(title);
  //TH2F *h1 = new TH2F("h1", "Histogram with Gaussian random distribution", fgNumOfPtBins, ptBinsArray, fgnumOfEtaBins, etaBinsArray);

  TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),fgNumOfPtBins, fgPtAxis, fgnumOfEtaBins, fgEtaAxis);
  histo->SetYTitle("#eta");
  histo->SetXTitle("p_{T}");
  histo->SetZTitle("E_{T}");
  histo->Sumw2();
  fhistoList->Add(histo);
  delete histoname;
  delete histotitle;
    
}

void AliAnalysisHadEt::CreateHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh)
{     //creates a 1d histogram of the given dimensions and adds it to the list of histograms to be saved
  TString *histoname   = new TString();
  TString *histotitle   = new TString();
  
  //cout<<"creating "<<name<<endl;

  histoname->Append(name);
  histotitle->Append(title);
  // printf("%s \n ",histoname->Data());
  TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
  histo->SetYTitle(ytitle);
  histo->SetXTitle(xtitle);
  histo->Sumw2();
  fhistoList->Add(histo);
  delete histoname;
  delete histotitle;
    
}
void AliAnalysisHadEt::CreateIntHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh)
{     //creates a 1d integer histogram and adds it to the list of histograms to be saved
  TString *histoname   = new TString();
  TString *histotitle   = new TString();
  
  //cout<<"creating "<<name<<endl;

  histoname->Append(name);
  histotitle->Append(title);
  // printf("%s \n ",histoname->Data());
  TH1I *histo = new TH1I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
  histo->SetYTitle(ytitle);
  histo->SetXTitle(xtitle);
  histo->Sumw2();
  fhistoList->Add(histo);
  delete histoname;
  delete histotitle;
    
}
void AliAnalysisHadEt::CreateHisto2D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh,Int_t ybins,Float_t ylow,Float_t yhigh)
{     //creates a 2d histogram and adds it to the list of histograms to be saved
  TString *histoname   = new TString();
  TString *histotitle   = new TString();
  
  //cout<<"creating "<<name<<endl;

  histoname->Append(name);
  histotitle->Append(title);
  // printf("%s \n ",histoname->Data());
  TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
  histo->SetYTitle(ytitle);
  histo->SetXTitle(xtitle);
  histo->Sumw2();
  fhistoList->Add(histo);
  delete histoname;
  delete histotitle;
    
}
void AliAnalysisHadEt::CreateIntHisto2D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh,Int_t ybins,Int_t ylow,Int_t yhigh)
{     //creates a 2-d integer histogram and adds it to the list of histograms to be saved
  TString *histoname   = new TString();
  TString *histotitle   = new TString();
  
  //cout<<"creating "<<name<<endl;

  histoname->Append(name);
  histotitle->Append(title);
  // printf("%s \n ",histoname->Data());
  TH2I *histo = new TH2I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
  histo->SetYTitle(ytitle);
  histo->SetXTitle(xtitle);
  histo->Sumw2();
  fhistoList->Add(histo);
  delete histoname;
  delete histotitle;
    
}

void AliAnalysisHadEt::CreateEtaHisto1D(TString name, TString title)
{     //creates 1d histogram in eta and adds it to the list of histograms to be saved
  TString *histoname   = new TString();
  TString *histotitle   = new TString();
  

  histoname->Append(name);
  histotitle->Append(title);
  TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),fgnumOfEtaBins, fgEtaAxis);
  histo->SetYTitle("E_{T}");
  histo->SetXTitle("#eta");
  histo->Sumw2();
  fhistoList->Add(histo);
  delete histoname;
  delete histotitle;
    
}
void AliAnalysisHadEt::FillHisto1D(TString histname, Float_t x, Float_t weight)
{//fills a 1d histogram with the name histoname with the value x and the weight "weight"
  TH1F     *histo; 
  TString  *name   = new TString();

  name->Append(histname);       
  histo = (TH1F *)fhistoList->FindObject(name->Data()); 
  if(histo){
    histo->Fill((Double_t)x, weight);
  }
  else{cerr<<"CorrelationMaker::FillHisto1D: no histogram "<< name->Data()<<endl;}
  delete name;
}
void AliAnalysisHadEt::FillHisto2D(TString histname, Float_t x, Float_t y, Float_t weight)
{//fills a 2d histogram with the name histoname with the value x and the weight "weight"
  TH2F     *histo; 
  TString  *name   = new TString();
  
  name->Append(histname);       
  histo = (TH2F *)fhistoList->FindObject(name->Data()); 
  if(histo){
    histo->Fill((Double_t)x,(Double_t)y, weight);
  }
  else{cerr<<"CorrelationMaker::FillHisto2D: no histogram "<< name->Data()<<endl;}
  delete name;
}


Float_t AliAnalysisHadEt::Et(TParticle *part, float mass){//function to calculate et in the same way as it would be calculated in a calorimeter
  if(mass+1000<0.01){//if no mass given return default.  The default argument is -1000
    if(TMath::Abs(part->GetPDG(0)->PdgCode())==2212 || TMath::Abs(part->GetPDG(0)->PdgCode())==2112){
      if(part->GetPDG(0)->PdgCode()==-2212 || part->GetPDG(0)->PdgCode()==-2112){//antiproton or antineutron
	//for antinucleons we specifically want to return the kinetic energy plus twice the rest mass
	return (part->Energy()+part->GetMass())*TMath::Sin(part->Theta());
      }
      if(part->GetPDG(0)->PdgCode()==2212 || part->GetPDG(0)->PdgCode()==2112){//antiproton or antineutron
	//for nucleons we specifically want to return the kinetic energy only
	return (part->Energy()-part->GetMass())*TMath::Sin(part->Theta());
      }
    }
    else{//otherwise go to the default
      return part->Energy()*TMath::Sin(part->Theta());
    }
  }
  else{//otherwise use the mass that was given
    return (TMath::Sqrt(TMath::Power(part->P(),2.0)+TMath::Power(mass,2.0)))*TMath::Sin(part->Theta());
  }
  return 0.0;
}
Float_t AliAnalysisHadEt::Et(Float_t p, Float_t theta, Int_t pid, Short_t charge){//function to calculate et in the same way as it would be calculated in a calorimeter
  if(pid==fPiPlusCode || pid==fPiMinusCode){//Nothing special for pions
    return TMath::Sqrt(p*p + fPionMass*fPionMass) * TMath::Sin(theta);
  }
  if(pid==fKPlusCode || pid==fKMinusCode){//Nothing special for kaons
    return TMath::Sqrt(p*p + fKaonMass*fKaonMass) * TMath::Sin(theta);
  }
  if(pid==fEPlusCode || pid==fEMinusCode){//Nothing special for electrons
    return TMath::Sqrt(p*p + fElectronMass*fElectronMass) * TMath::Sin(theta);
  }
  if(pid==fProtonCode || pid==fAntiProtonCode){//But for protons we must be careful...
    if(charge<0.0){//antiprotns: kinetic energy plus twice the rest mass
      return (TMath::Sqrt(p*p + fProtonMass*fProtonMass) + fProtonMass) * TMath::Sin(theta);
    }
    if(charge>0.0){//antiprotns: kinetic energy only
      return (TMath::Sqrt(p*p + fProtonMass*fProtonMass) - fProtonMass) * TMath::Sin(theta);
    }
  }
  cerr<<"Uh-oh!  Et not set properly!"<<endl;
  return 0.0;
}
Commit	Line	Data
cf6522d1	1	//Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
641e1e0c	2	//University of Tennessee at Knoxville
cf6522d1	3	//
	4	// This class is designed for the analysis of the hadronic component of
	5	// transverse energy. It is used by AliAnalysisTaskHadEt.
	6	// This gets information about the hadronic component of the transverse energy
	7	// from tracks reconstructed in an event
	8	// it has daughters, AliAnalysisHadEtMonteCarlo and
	9	// AliAnalysisHadEtReconstructed which loop over either Monte Carlo data or
	10	// real data to get Et
3ce6b879	11
641e1e0c	12	#include "AliAnalysisHadEt.h"
	13	#include "TMath.h"
	14	#include "TList.h"
	15	#include "TH1F.h"
	16	#include "TH2F.h"
	17	#include <iostream>
	18	#include "AliAnalysisEtCuts.h"
020abff0	19	#include "AliMCEvent.h"
641e1e0c	20	#include "AliVEvent.h"
020abff0	21	#include "AliStack.h"
	22	#include "AliESDtrackCuts.h"
	23	#include "TDatabasePDG.h"
	24	#include "TParticle.h"
	25	#include "Rtypes.h"
641e1e0c	26
	27	using namespace std;
	28
	29	ClassImp(AliAnalysisHadEt);
	30
	31
020abff0	32	Int_t AliAnalysisHadEt::fgnumOfEtaBins = 46;
	33	Float_t AliAnalysisHadEt::fgEtaAxis[47]={-0.78, -0.74, -0.7, -0.66, -0.62, -0.58, -0.54, -0.5, -0.46, -0.42, -0.38, -0.34, -0.3, -0.26, -0.22, -0.18, -0.14, -0.12, -0.1, -0.08, -0.06, -0.04, -0.02, -0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.18, 0.22, 0.26, 0.3, 0.34, 0.38, 0.42, 0.46, 0.5, 0.54, 0.58, 0.62, 0.66, 0.7, 0.74, 0.78};
	34	Int_t AliAnalysisHadEt::fgNumOfPtBins = 111;
	35	Float_t AliAnalysisHadEt::fgPtAxis[117]=
641e1e0c	36	{0.0,0.01,0.02,0.03,0.04, 0.05, 0.06,0.07,0.08,0.09, 0.10,0.11, .12,0.13, .14,0.15, .16,0.17, .18,0.19,
	37	0.2, .22, .24, .26, .28, 0.30, 0.32, .34, .36, .38, 0.40, .42, .44, .46, .48,
	38	0.5, .52, .54, .56, .58, 0.60, 0.62, .64, .66, .68, 0.70, .72, .74, .76, .78,
	39	.80, .82, .84, .86, .88, 0.90, 0.92, .94, .96, .98, 1.00,1.05, 1.1,1.15, 1.2,
	40	1.25, 1.3,1.35,1.40,1.45, 1.50, 1.55, 1.6,1.65, 1.7, 1.75, 1.8,1.85, 1.9,1.95,
	41	2.0, 2.2, 2.4, 2.6, 2.8, 3.00, 3.20, 3.4, 3.6, 3.8, 4.00, 4.2, 4.4, 4.6, 4.8,
	42	5.0, 5.5, 6.0, 6.5, 7.0, 7.50, 8.00, 8.5, 9.0, 9.5, 10.0,12.0,14.0,16.0,18.0,
	43	20.0,25.0,30.0,35.0,40.0, 45.0, 50.0};
	44
d6214a64	45	Float_t AliAnalysisHadEt::fgPtTPCCutOff = 0.15;
	46	Float_t AliAnalysisHadEt::fgPtITSCutOff = 0.10;
	47
641e1e0c	48	AliAnalysisHadEt::AliAnalysisHadEt() :
641e1e0c	49	fHistogramNameSuffix("")
4998becf	50	,fCuts(0)
641e1e0c	51	,fPdgDB(0)
020abff0	52	,fPiPlusCode(0)
	53	,fPiMinusCode(0)
	54	,fKPlusCode(0)
	55	,fKMinusCode(0)
	56	,fProtonCode(0)
	57	,fAntiProtonCode(0)
	58	,fLambdaCode(0)
	59	,fAntiLambdaCode(0)
	60	,fK0SCode(0)
	61	,fOmegaCode(0)
	62	,fAntiOmegaCode(0)
	63	,fXi0Code(0)
	64	,fAntiXi0Code(0)
	65	,fXiCode(0)
	66	,fAntiXiCode(0)
	67	,fSigmaCode(0)
	68	,fAntiSigmaCode(0)
	69	,fK0LCode(0)
	70	,fNeutronCode(0)
	71	,fAntiNeutronCode(0)
	72	,fEPlusCode(0)
	73	,fEMinusCode(0)
f427cbed	74	,fGammaCode(0)
	75	,fPi0Code(0)
	76	,fEtaCode(0)
b64de20c	77	,fOmega0Code(0)
020abff0	78	,fPionMass(0)
3ce6b879	79	,fKaonMass(0)
	80	,fProtonMass(0)
	81	,fElectronMass(0)
641e1e0c	82	,fSumEt(0)
	83	,fSumEtAcc(0)
	84	,fTotEt(0)
	85	,fTotEtAcc(0)
	86	,fTotNeutralEt(0)
	87	,fTotNeutralEtAcc(0)
	88	,fTotChargedEt(0)
	89	,fTotChargedEtAcc(0)
	90	,fMultiplicity(0)
	91	,fChargedMultiplicity(0)
	92	,fNeutralMultiplicity(0)
cf6522d1	93	,fEsdtrackCutsITSTPC(0)
	94	,fEsdtrackCutsTPC(0)
	95	,fEsdtrackCutsITS(0)
020abff0	96	,fhistoList(0)
020abff0	97	{//default constructor
641e1e0c	98
	99	}
	100
	101	AliAnalysisHadEt::~AliAnalysisHadEt()
020abff0	102	{//destructor
641e1e0c	103
	104	}
	105
	106	Int_t AliAnalysisHadEt::AnalyseEvent(AliVEvent *event)
020abff0	107	{ //this line is basically here to eliminate a compiler warning that event is not used. Making it a virtual function did not work with the plugin.
0065531a	108	cout << "This event has " << event->GetNumberOfTracks() << " tracks" << endl;
0065531a	109	ResetEventValues();
641e1e0c	110	return 0;
	111	}
	112
	113	void AliAnalysisHadEt::FillOutputList()
020abff0	114	{//fill the output histogram list with histograms in all AliAnalysisHadEt's
641e1e0c	115	}
	116
	117	void AliAnalysisHadEt::Init()
0065531a	118	{// clear variables, set up cuts and PDG info
0065531a	119	ResetEventValues();
3ce6b879	120
641e1e0c	121	}
	122
	123	void AliAnalysisHadEt::CreateHistograms()
020abff0	124	{//creates histograms included in all AliAnalysisHadEt's
641e1e0c	125	}
	126
	127	void AliAnalysisHadEt::FillHistograms()
020abff0	128	{//Fills histograms filled for all AliAnalysisHadEt's
641e1e0c	129	}
	130
	131	void AliAnalysisHadEt::ResetEventValues()
020abff0	132	{//Resets event values of et to zero
0065531a	133	fTotEt = 0;
	134	fTotEtAcc = 0;
	135	fTotNeutralEt = 0;
	136	fTotNeutralEtAcc = 0;
	137	fTotChargedEt = 0;
	138	fTotChargedEtAcc = 0;
	139	fMultiplicity = 0;
	140	fChargedMultiplicity = 0;
	141	fNeutralMultiplicity = 0;
	142
	143	if (!fCuts \|\| !fPdgDB \|\| fPiPlusCode==0) { // some Init's needed
	144	cout << __FILE__ << ":" << __LINE__ << " : Init " << endl;
	145	if (!fCuts) {
	146	cout << " setting up Cuts " << endl;
	147	fCuts = new AliAnalysisEtCuts();
	148	}
	149	if(!fPdgDB) {
	150	cout << " setting up PdgDB " << endl;
	151	fPdgDB = new TDatabasePDG();
	152	}
	153
	154	if (fPiPlusCode==0) {
	155	SetParticleCodes();
	156	}
	157	}
641e1e0c	158	}
4998becf	159
	160	void AliAnalysisHadEt::SetParticleCodes()
	161	{ //the codes are defined in $ROOTSYS/etc/pdg_table.txt
f427cbed	162	fPionMass = fPdgDB->GetParticle("pi+")->Mass();
3ce6b879	163	fKaonMass = fPdgDB->GetParticle("K+")->Mass();
	164	fProtonMass = fPdgDB->GetParticle("proton")->Mass();
	165	fElectronMass = fPdgDB->GetParticle("e+")->Mass();
f427cbed	166	fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
	167	fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
	168	fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
	169	fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
	170	fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
	171	fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
	172	fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
	173	fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
	174	fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
	175	fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
	176	fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
	177	fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
	178	fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
	179	fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
	180	fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
	181	fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
	182	fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
	183	fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
	184	fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
	185	fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
	186	fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
	187	fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
	188	fGammaCode = fPdgDB->GetParticle("gamma")->PdgCode();
	189	fPi0Code = fPdgDB->GetParticle("pi0")->PdgCode();
	190	fEtaCode = fPdgDB->GetParticle("eta")->PdgCode();
b64de20c	191	fOmega0Code = fPdgDB->GetParticle("omega")->PdgCode();
f427cbed	192	cout << "Resetting Codes: Pion " << fPiPlusCode
	193	<< "," << fPiMinusCode
	194	<< " Kaon " << fKPlusCode
	195	<< "," << fKMinusCode
	196	<< " Gamma "<<fGammaCode
	197	<< " Pi0 "<<fPi0Code
	198	<< " Eta "<<fEtaCode
b64de20c	199	<< " omega "<<fOmega0Code
f427cbed	200	<< endl;
4998becf	201	}
4998becf	202
641e1e0c	203	void AliAnalysisHadEt::CreateEtaPtHisto2D(TString name, TString title)
020abff0	204	{ //creates a 2-d histogram in eta and phi and adds it to the list of histograms to be saved
641e1e0c	205	TString *histoname = new TString();
	206	TString *histotitle = new TString();
	207
	208	histoname->Append(name);
	209	histotitle->Append(title);
020abff0	210	//TH2F *h1 = new TH2F("h1", "Histogram with Gaussian random distribution", fgNumOfPtBins, ptBinsArray, fgnumOfEtaBins, etaBinsArray);
641e1e0c	211
020abff0	212	TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),fgNumOfPtBins, fgPtAxis, fgnumOfEtaBins, fgEtaAxis);
641e1e0c	213	histo->SetYTitle("#eta");
	214	histo->SetXTitle("p_{T}");
	215	histo->SetZTitle("E_{T}");
	216	histo->Sumw2();
020abff0	217	fhistoList->Add(histo);
641e1e0c	218	delete histoname;
	219	delete histotitle;
	220
	221	}
	222
	223	void AliAnalysisHadEt::CreateHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh)
020abff0	224	{ //creates a 1d histogram of the given dimensions and adds it to the list of histograms to be saved
641e1e0c	225	TString *histoname = new TString();
	226	TString *histotitle = new TString();
	227
	228	//cout<<"creating "<<name<<endl;
	229
	230	histoname->Append(name);
	231	histotitle->Append(title);
	232	// printf("%s \n ",histoname->Data());
	233	TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
	234	histo->SetYTitle(ytitle);
	235	histo->SetXTitle(xtitle);
	236	histo->Sumw2();
020abff0	237	fhistoList->Add(histo);
641e1e0c	238	delete histoname;
	239	delete histotitle;
	240
	241	}
	242	void AliAnalysisHadEt::CreateIntHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh)
020abff0	243	{ //creates a 1d integer histogram and adds it to the list of histograms to be saved
641e1e0c	244	TString *histoname = new TString();
	245	TString *histotitle = new TString();
	246
	247	//cout<<"creating "<<name<<endl;
	248
	249	histoname->Append(name);
	250	histotitle->Append(title);
	251	// printf("%s \n ",histoname->Data());
	252	TH1I *histo = new TH1I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
	253	histo->SetYTitle(ytitle);
	254	histo->SetXTitle(xtitle);
	255	histo->Sumw2();
020abff0	256	fhistoList->Add(histo);
641e1e0c	257	delete histoname;
	258	delete histotitle;
	259
	260	}
	261	void AliAnalysisHadEt::CreateHisto2D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh,Int_t ybins,Float_t ylow,Float_t yhigh)
020abff0	262	{ //creates a 2d histogram and adds it to the list of histograms to be saved
641e1e0c	263	TString *histoname = new TString();
	264	TString *histotitle = new TString();
	265
	266	//cout<<"creating "<<name<<endl;
	267
	268	histoname->Append(name);
	269	histotitle->Append(title);
	270	// printf("%s \n ",histoname->Data());
	271	TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
	272	histo->SetYTitle(ytitle);
	273	histo->SetXTitle(xtitle);
	274	histo->Sumw2();
020abff0	275	fhistoList->Add(histo);
641e1e0c	276	delete histoname;
	277	delete histotitle;
	278
	279	}
	280	void AliAnalysisHadEt::CreateIntHisto2D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh,Int_t ybins,Int_t ylow,Int_t yhigh)
020abff0	281	{ //creates a 2-d integer histogram and adds it to the list of histograms to be saved
641e1e0c	282	TString *histoname = new TString();
	283	TString *histotitle = new TString();
	284
	285	//cout<<"creating "<<name<<endl;
	286
	287	histoname->Append(name);
	288	histotitle->Append(title);
	289	// printf("%s \n ",histoname->Data());
	290	TH2I *histo = new TH2I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
	291	histo->SetYTitle(ytitle);
	292	histo->SetXTitle(xtitle);
	293	histo->Sumw2();
020abff0	294	fhistoList->Add(histo);
641e1e0c	295	delete histoname;
	296	delete histotitle;
	297
	298	}
	299
	300	void AliAnalysisHadEt::CreateEtaHisto1D(TString name, TString title)
020abff0	301	{ //creates 1d histogram in eta and adds it to the list of histograms to be saved
641e1e0c	302	TString *histoname = new TString();
	303	TString *histotitle = new TString();
	304
	305
	306	histoname->Append(name);
	307	histotitle->Append(title);
020abff0	308	TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),fgnumOfEtaBins, fgEtaAxis);
641e1e0c	309	histo->SetYTitle("E_{T}");
	310	histo->SetXTitle("#eta");
	311	histo->Sumw2();
020abff0	312	fhistoList->Add(histo);
641e1e0c	313	delete histoname;
	314	delete histotitle;
	315
	316	}
	317	void AliAnalysisHadEt::FillHisto1D(TString histname, Float_t x, Float_t weight)
020abff0	318	{//fills a 1d histogram with the name histoname with the value x and the weight "weight"
641e1e0c	319	TH1F *histo;
	320	TString *name = new TString();
	321
	322	name->Append(histname);
020abff0	323	histo = (TH1F *)fhistoList->FindObject(name->Data());
641e1e0c	324	if(histo){
	325	histo->Fill((Double_t)x, weight);
	326	}
	327	else{cerr<<"CorrelationMaker::FillHisto1D: no histogram "<< name->Data()<<endl;}
	328	delete name;
	329	}
	330	void AliAnalysisHadEt::FillHisto2D(TString histname, Float_t x, Float_t y, Float_t weight)
020abff0	331	{//fills a 2d histogram with the name histoname with the value x and the weight "weight"
641e1e0c	332	TH2F *histo;
	333	TString *name = new TString();
	334
	335	name->Append(histname);
020abff0	336	histo = (TH2F *)fhistoList->FindObject(name->Data());
641e1e0c	337	if(histo){
	338	histo->Fill((Double_t)x,(Double_t)y, weight);
	339	}
	340	else{cerr<<"CorrelationMaker::FillHisto2D: no histogram "<< name->Data()<<endl;}
	341	delete name;
	342	}
	343
	344
020abff0	345	Float_t AliAnalysisHadEt::Et(TParticle *part, float mass){//function to calculate et in the same way as it would be calculated in a calorimeter
020abff0	346	if(mass+1000<0.01){//if no mass given return default. The default argument is -1000
641e1e0c	347	if(TMath::Abs(part->GetPDG(0)->PdgCode())==2212 \|\| TMath::Abs(part->GetPDG(0)->PdgCode())==2112){
	348	if(part->GetPDG(0)->PdgCode()==-2212 \|\| part->GetPDG(0)->PdgCode()==-2112){//antiproton or antineutron
	349	//for antinucleons we specifically want to return the kinetic energy plus twice the rest mass
	350	return (part->Energy()+part->GetMass())*TMath::Sin(part->Theta());
	351	}
	352	if(part->GetPDG(0)->PdgCode()==2212 \|\| part->GetPDG(0)->PdgCode()==2112){//antiproton or antineutron
	353	//for nucleons we specifically want to return the kinetic energy only
	354	return (part->Energy()-part->GetMass())*TMath::Sin(part->Theta());
	355	}
	356	}
	357	else{//otherwise go to the default
	358	return part->Energy()*TMath::Sin(part->Theta());
	359	}
	360	}
	361	else{//otherwise use the mass that was given
	362	return (TMath::Sqrt(TMath::Power(part->P(),2.0)+TMath::Power(mass,2.0)))*TMath::Sin(part->Theta());
	363	}
	364	return 0.0;
	365	}
3ce6b879	366	Float_t AliAnalysisHadEt::Et(Float_t p, Float_t theta, Int_t pid, Short_t charge){//function to calculate et in the same way as it would be calculated in a calorimeter
	367	if(pid==fPiPlusCode \|\| pid==fPiMinusCode){//Nothing special for pions
	368	return TMath::Sqrt(pp + fPionMassfPionMass) * TMath::Sin(theta);
	369	}
	370	if(pid==fKPlusCode \|\| pid==fKMinusCode){//Nothing special for kaons
	371	return TMath::Sqrt(pp + fKaonMassfKaonMass) * TMath::Sin(theta);
	372	}
	373	if(pid==fEPlusCode \|\| pid==fEMinusCode){//Nothing special for electrons
	374	return TMath::Sqrt(pp + fElectronMassfElectronMass) * TMath::Sin(theta);
	375	}
	376	if(pid==fProtonCode \|\| pid==fAntiProtonCode){//But for protons we must be careful...
	377	if(charge<0.0){//antiprotns: kinetic energy plus twice the rest mass
	378	return (TMath::Sqrt(pp + fProtonMassfProtonMass) + fProtonMass) * TMath::Sin(theta);
	379	}
	380	if(charge>0.0){//antiprotns: kinetic energy only
	381	return (TMath::Sqrt(pp + fProtonMassfProtonMass) - fProtonMass) * TMath::Sin(theta);
	382	}
	383	}
	384	cerr<<"Uh-oh! Et not set properly!"<<endl;
	385	return 0.0;
	386	}