[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}; 

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
	45	AliAnalysisHadEt::AliAnalysisHadEt() :
	46	fHistogramNameSuffix("")
4998becf	47	,fCuts(0)
641e1e0c	48	,fPdgDB(0)
020abff0	49	,fPiPlusCode(0)
	50	,fPiMinusCode(0)
	51	,fKPlusCode(0)
	52	,fKMinusCode(0)
	53	,fProtonCode(0)
	54	,fAntiProtonCode(0)
	55	,fLambdaCode(0)
	56	,fAntiLambdaCode(0)
	57	,fK0SCode(0)
	58	,fOmegaCode(0)
	59	,fAntiOmegaCode(0)
	60	,fXi0Code(0)
	61	,fAntiXi0Code(0)
	62	,fXiCode(0)
	63	,fAntiXiCode(0)
	64	,fSigmaCode(0)
	65	,fAntiSigmaCode(0)
	66	,fK0LCode(0)
	67	,fNeutronCode(0)
	68	,fAntiNeutronCode(0)
	69	,fEPlusCode(0)
	70	,fEMinusCode(0)
f427cbed	71	,fGammaCode(0)
	72	,fPi0Code(0)
	73	,fEtaCode(0)
b64de20c	74	,fOmega0Code(0)
020abff0	75	,fPionMass(0)
3ce6b879	76	,fKaonMass(0)
	77	,fProtonMass(0)
	78	,fElectronMass(0)
641e1e0c	79	,fSumEt(0)
	80	,fSumEtAcc(0)
	81	,fTotEt(0)
	82	,fTotEtAcc(0)
	83	,fTotNeutralEt(0)
	84	,fTotNeutralEtAcc(0)
	85	,fTotChargedEt(0)
	86	,fTotChargedEtAcc(0)
	87	,fMultiplicity(0)
	88	,fChargedMultiplicity(0)
	89	,fNeutralMultiplicity(0)
cf6522d1	90	,fEsdtrackCutsITSTPC(0)
	91	,fEsdtrackCutsTPC(0)
	92	,fEsdtrackCutsITS(0)
020abff0	93	,fhistoList(0)
020abff0	94	{//default constructor
641e1e0c	95
	96	}
	97
	98	AliAnalysisHadEt::~AliAnalysisHadEt()
020abff0	99	{//destructor
641e1e0c	100
	101	}
	102
	103	Int_t AliAnalysisHadEt::AnalyseEvent(AliVEvent *event)
020abff0	104	{ //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	105	cout << "This event has " << event->GetNumberOfTracks() << " tracks" << endl;
0065531a	106	ResetEventValues();
641e1e0c	107	return 0;
	108	}
	109
	110	void AliAnalysisHadEt::FillOutputList()
020abff0	111	{//fill the output histogram list with histograms in all AliAnalysisHadEt's
641e1e0c	112	}
	113
	114	void AliAnalysisHadEt::Init()
0065531a	115	{// clear variables, set up cuts and PDG info
0065531a	116	ResetEventValues();
3ce6b879	117
641e1e0c	118	}
	119
	120	void AliAnalysisHadEt::CreateHistograms()
020abff0	121	{//creates histograms included in all AliAnalysisHadEt's
641e1e0c	122	}
	123
	124	void AliAnalysisHadEt::FillHistograms()
020abff0	125	{//Fills histograms filled for all AliAnalysisHadEt's
641e1e0c	126	}
	127
	128	void AliAnalysisHadEt::ResetEventValues()
020abff0	129	{//Resets event values of et to zero
0065531a	130	fTotEt = 0;
	131	fTotEtAcc = 0;
	132	fTotNeutralEt = 0;
	133	fTotNeutralEtAcc = 0;
	134	fTotChargedEt = 0;
	135	fTotChargedEtAcc = 0;
	136	fMultiplicity = 0;
	137	fChargedMultiplicity = 0;
	138	fNeutralMultiplicity = 0;
	139
	140	if (!fCuts \|\| !fPdgDB \|\| fPiPlusCode==0) { // some Init's needed
	141	cout << __FILE__ << ":" << __LINE__ << " : Init " << endl;
	142	if (!fCuts) {
	143	cout << " setting up Cuts " << endl;
	144	fCuts = new AliAnalysisEtCuts();
	145	}
	146	if(!fPdgDB) {
	147	cout << " setting up PdgDB " << endl;
	148	fPdgDB = new TDatabasePDG();
	149	}
	150
	151	if (fPiPlusCode==0) {
	152	SetParticleCodes();
	153	}
	154	}
641e1e0c	155	}
4998becf	156
	157	void AliAnalysisHadEt::SetParticleCodes()
	158	{ //the codes are defined in $ROOTSYS/etc/pdg_table.txt
f427cbed	159	fPionMass = fPdgDB->GetParticle("pi+")->Mass();
3ce6b879	160	fKaonMass = fPdgDB->GetParticle("K+")->Mass();
	161	fProtonMass = fPdgDB->GetParticle("proton")->Mass();
	162	fElectronMass = fPdgDB->GetParticle("e+")->Mass();
f427cbed	163	fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
	164	fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
	165	fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
	166	fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
	167	fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
	168	fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
	169	fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
	170	fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
	171	fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
	172	fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
	173	fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
	174	fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
	175	fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
	176	fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
	177	fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
	178	fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
	179	fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
	180	fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
	181	fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
	182	fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
	183	fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
	184	fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
	185	fGammaCode = fPdgDB->GetParticle("gamma")->PdgCode();
	186	fPi0Code = fPdgDB->GetParticle("pi0")->PdgCode();
	187	fEtaCode = fPdgDB->GetParticle("eta")->PdgCode();
b64de20c	188	fOmega0Code = fPdgDB->GetParticle("omega")->PdgCode();
f427cbed	189	cout << "Resetting Codes: Pion " << fPiPlusCode
	190	<< "," << fPiMinusCode
	191	<< " Kaon " << fKPlusCode
	192	<< "," << fKMinusCode
	193	<< " Gamma "<<fGammaCode
	194	<< " Pi0 "<<fPi0Code
	195	<< " Eta "<<fEtaCode
b64de20c	196	<< " omega "<<fOmega0Code
f427cbed	197	<< endl;
4998becf	198	}
4998becf	199
641e1e0c	200	void AliAnalysisHadEt::CreateEtaPtHisto2D(TString name, TString title)
020abff0	201	{ //creates a 2-d histogram in eta and phi and adds it to the list of histograms to be saved
641e1e0c	202	TString *histoname = new TString();
	203	TString *histotitle = new TString();
	204
	205	histoname->Append(name);
	206	histotitle->Append(title);
020abff0	207	//TH2F *h1 = new TH2F("h1", "Histogram with Gaussian random distribution", fgNumOfPtBins, ptBinsArray, fgnumOfEtaBins, etaBinsArray);
641e1e0c	208
020abff0	209	TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),fgNumOfPtBins, fgPtAxis, fgnumOfEtaBins, fgEtaAxis);
641e1e0c	210	histo->SetYTitle("#eta");
	211	histo->SetXTitle("p_{T}");
	212	histo->SetZTitle("E_{T}");
	213	histo->Sumw2();
020abff0	214	fhistoList->Add(histo);
641e1e0c	215	delete histoname;
	216	delete histotitle;
	217
	218	}
	219
	220	void AliAnalysisHadEt::CreateHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh)
020abff0	221	{ //creates a 1d histogram of the given dimensions and adds it to the list of histograms to be saved
641e1e0c	222	TString *histoname = new TString();
	223	TString *histotitle = new TString();
	224
	225	//cout<<"creating "<<name<<endl;
	226
	227	histoname->Append(name);
	228	histotitle->Append(title);
	229	// printf("%s \n ",histoname->Data());
	230	TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
	231	histo->SetYTitle(ytitle);
	232	histo->SetXTitle(xtitle);
	233	histo->Sumw2();
020abff0	234	fhistoList->Add(histo);
641e1e0c	235	delete histoname;
	236	delete histotitle;
	237
	238	}
	239	void AliAnalysisHadEt::CreateIntHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh)
020abff0	240	{ //creates a 1d integer histogram and adds it to the list of histograms to be saved
641e1e0c	241	TString *histoname = new TString();
	242	TString *histotitle = new TString();
	243
	244	//cout<<"creating "<<name<<endl;
	245
	246	histoname->Append(name);
	247	histotitle->Append(title);
	248	// printf("%s \n ",histoname->Data());
	249	TH1I *histo = new TH1I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
	250	histo->SetYTitle(ytitle);
	251	histo->SetXTitle(xtitle);
	252	histo->Sumw2();
020abff0	253	fhistoList->Add(histo);
641e1e0c	254	delete histoname;
	255	delete histotitle;
	256
	257	}
	258	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	259	{ //creates a 2d histogram and adds it to the list of histograms to be saved
641e1e0c	260	TString *histoname = new TString();
	261	TString *histotitle = new TString();
	262
	263	//cout<<"creating "<<name<<endl;
	264
	265	histoname->Append(name);
	266	histotitle->Append(title);
	267	// printf("%s \n ",histoname->Data());
	268	TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
	269	histo->SetYTitle(ytitle);
	270	histo->SetXTitle(xtitle);
	271	histo->Sumw2();
020abff0	272	fhistoList->Add(histo);
641e1e0c	273	delete histoname;
	274	delete histotitle;
	275
	276	}
	277	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	278	{ //creates a 2-d integer histogram and adds it to the list of histograms to be saved
641e1e0c	279	TString *histoname = new TString();
	280	TString *histotitle = new TString();
	281
	282	//cout<<"creating "<<name<<endl;
	283
	284	histoname->Append(name);
	285	histotitle->Append(title);
	286	// printf("%s \n ",histoname->Data());
	287	TH2I *histo = new TH2I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
	288	histo->SetYTitle(ytitle);
	289	histo->SetXTitle(xtitle);
	290	histo->Sumw2();
020abff0	291	fhistoList->Add(histo);
641e1e0c	292	delete histoname;
	293	delete histotitle;
	294
	295	}
	296
	297	void AliAnalysisHadEt::CreateEtaHisto1D(TString name, TString title)
020abff0	298	{ //creates 1d histogram in eta and adds it to the list of histograms to be saved
641e1e0c	299	TString *histoname = new TString();
	300	TString *histotitle = new TString();
	301
	302
	303	histoname->Append(name);
	304	histotitle->Append(title);
020abff0	305	TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),fgnumOfEtaBins, fgEtaAxis);
641e1e0c	306	histo->SetYTitle("E_{T}");
	307	histo->SetXTitle("#eta");
	308	histo->Sumw2();
020abff0	309	fhistoList->Add(histo);
641e1e0c	310	delete histoname;
	311	delete histotitle;
	312
	313	}
	314	void AliAnalysisHadEt::FillHisto1D(TString histname, Float_t x, Float_t weight)
020abff0	315	{//fills a 1d histogram with the name histoname with the value x and the weight "weight"
641e1e0c	316	TH1F *histo;
	317	TString *name = new TString();
	318
	319	name->Append(histname);
020abff0	320	histo = (TH1F *)fhistoList->FindObject(name->Data());
641e1e0c	321	if(histo){
	322	histo->Fill((Double_t)x, weight);
	323	}
	324	else{cerr<<"CorrelationMaker::FillHisto1D: no histogram "<< name->Data()<<endl;}
	325	delete name;
	326	}
	327	void AliAnalysisHadEt::FillHisto2D(TString histname, Float_t x, Float_t y, Float_t weight)
020abff0	328	{//fills a 2d histogram with the name histoname with the value x and the weight "weight"
641e1e0c	329	TH2F *histo;
	330	TString *name = new TString();
	331
	332	name->Append(histname);
020abff0	333	histo = (TH2F *)fhistoList->FindObject(name->Data());
641e1e0c	334	if(histo){
	335	histo->Fill((Double_t)x,(Double_t)y, weight);
	336	}
	337	else{cerr<<"CorrelationMaker::FillHisto2D: no histogram "<< name->Data()<<endl;}
	338	delete name;
	339	}
	340
	341
020abff0	342	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	343	if(mass+1000<0.01){//if no mass given return default. The default argument is -1000
641e1e0c	344	if(TMath::Abs(part->GetPDG(0)->PdgCode())==2212 \|\| TMath::Abs(part->GetPDG(0)->PdgCode())==2112){
	345	if(part->GetPDG(0)->PdgCode()==-2212 \|\| part->GetPDG(0)->PdgCode()==-2112){//antiproton or antineutron
	346	//for antinucleons we specifically want to return the kinetic energy plus twice the rest mass
	347	return (part->Energy()+part->GetMass())*TMath::Sin(part->Theta());
	348	}
	349	if(part->GetPDG(0)->PdgCode()==2212 \|\| part->GetPDG(0)->PdgCode()==2112){//antiproton or antineutron
	350	//for nucleons we specifically want to return the kinetic energy only
	351	return (part->Energy()-part->GetMass())*TMath::Sin(part->Theta());
	352	}
	353	}
	354	else{//otherwise go to the default
	355	return part->Energy()*TMath::Sin(part->Theta());
	356	}
	357	}
	358	else{//otherwise use the mass that was given
	359	return (TMath::Sqrt(TMath::Power(part->P(),2.0)+TMath::Power(mass,2.0)))*TMath::Sin(part->Theta());
	360	}
	361	return 0.0;
	362	}
3ce6b879	363	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
	364	if(pid==fPiPlusCode \|\| pid==fPiMinusCode){//Nothing special for pions
	365	return TMath::Sqrt(pp + fPionMassfPionMass) * TMath::Sin(theta);
	366	}
	367	if(pid==fKPlusCode \|\| pid==fKMinusCode){//Nothing special for kaons
	368	return TMath::Sqrt(pp + fKaonMassfKaonMass) * TMath::Sin(theta);
	369	}
	370	if(pid==fEPlusCode \|\| pid==fEMinusCode){//Nothing special for electrons
	371	return TMath::Sqrt(pp + fElectronMassfElectronMass) * TMath::Sin(theta);
	372	}
	373	if(pid==fProtonCode \|\| pid==fAntiProtonCode){//But for protons we must be careful...
	374	if(charge<0.0){//antiprotns: kinetic energy plus twice the rest mass
	375	return (TMath::Sqrt(pp + fProtonMassfProtonMass) + fProtonMass) * TMath::Sin(theta);
	376	}
	377	if(charge>0.0){//antiprotns: kinetic energy only
	378	return (TMath::Sqrt(pp + fProtonMassfProtonMass) - fProtonMass) * TMath::Sin(theta);
	379	}
	380	}
	381	cerr<<"Uh-oh! Et not set properly!"<<endl;
	382	return 0.0;
	383	}