[u/mrichter/AliRoot.git] / PWG0 / multPbPb / AliAnalysisMultPbTrackHistoManager.cxx

#include "AliAnalysisMultPbTrackHistoManager.h"
#include "AliLog.h"
#include "TH1.h"
#include "TH3D.h"
#include "TH1I.h"
#include "TROOT.h"
#include "TMCProcess.h"
#include "AliMCParticle.h"

#include <iostream>
#include "TH2D.h"

using namespace std;

ClassImp(AliAnalysisMultPbTrackHistoManager)

const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[]     = { "All Events", "After physics Selection", "After centrality selection", "With Vertex" };
const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim", 
									  "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[]     = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[]     = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake"};
const char * AliAnalysisMultPbTrackHistoManager::kSpeciesName[]     = { "pi+", "K+", "p", "l+",  "pi-", "K-", "barp", "l-", "Others"};


AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager() : AliHistoListWrapper(), fHNameSuffix(""){ 
  // standard ctor

}

AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const char * name, const char * title): AliHistoListWrapper(name,title), fHNameSuffix("")  {
  //named ctor
};

AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) : AliHistoListWrapper (obj) {
  // copy ctor
  AliError("Not Implemented");
};

AliAnalysisMultPbTrackHistoManager::~AliAnalysisMultPbTrackHistoManager() {
  // dtor

}

TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id, Int_t particle) {
  // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.

  TString name = kHistoPtEtaVzNames[id];

  if(particle >= 0) name += kSpeciesName[particle];

  TH3D * h = (TH3D*) GetHisto(name.Data());
  if (!h) {
    h = BookHistoPtEtaVz(name.Data(), Form("Pt Eta Vz distribution (%s)",kHistoPtEtaVzNames[id]));
  }

  return h;

}

TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
  // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.

  TH1D * h = (TH1D*) GetHisto(kHistoDCANames[id]);
  if (!h) {
    h = BookHistoDCA(kHistoDCANames[id], Form("Pt Eta Vz distribution (%s)",kHistoDCANames[id]));
  }

  return h;

}

TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMult(Histo_t id) {
  // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.

  TString name = kHistoPrefix[id];
  name += "Mult";
  TH1D * h = (TH1D*) GetHisto(name.Data());
  if (!h) {
    h = BookHistoMult(name.Data(), Form("Multiplicity distribution (%s)",kHistoPrefix[id]));
  }

  return h;

}

TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id, 
						       Float_t minEta, Float_t maxEta, 
						       Float_t minVz, Float_t maxVz, 
						       Bool_t scaleWidth) {
  // Returns a projection of the 3D histo pt/eta/vz.
  // WARNING: since that is a histo, the requested range will be discretized to the binning.
  // Always avoids under (over) flows
  // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)

  TH3D * h3D = GetHistoPtEtaVz(id);

  // Get range in terms of bin numners.  If the float range is
  // less than -11111 take the range from the first to the last bin (i.e. no
  // under/over-flows)
  // FIXME: UNDERFLOWS
  Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
  Int_t min2  = minVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;

  // Int_t max1 = maxEta  < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
  // Int_t max2  = maxVz  < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
  Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
  Int_t max2  = maxVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(maxVz -0.00001);


  TString hname = h3D->GetName();
  hname = hname +  "_pt_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);

  
  if (gROOT->FindObjectAny(hname.Data())){
    AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
    hname += "_2";
  }

  TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E");
  if(scaleWidth) h ->Scale(1.,"width");

  return h;

}


TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoPtVz (Histo_t id, 
							 Float_t minEta, Float_t maxEta, 
							 Bool_t scaleWidth) {
  // Returns a projection of the 3D histo pt/eta/vz.
  // WARNING: since that is a histo, the requested range will be discretized to the binning.
  // Always avoids under (over) flows
  // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)

  // FIXME: what do I do here for the scaling?

  TH3D * h3D = GetHistoPtEtaVz(id);

  // Get range in terms of bin numners.  If the float range is
  // less than -11111 take the range from the first to the last bin (i.e. no
  // under/over-flows)
  Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
  Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);


  TString hname = h3D->GetName();
  hname = hname +  "_ptvz_" + long (min1)  + "_" + long(max1);
  
  if (gROOT->FindObjectAny(hname.Data())){
    AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
    hname += "_2";
  }

  h3D->GetYaxis()->SetRange(min1,max1);
  
  TH2D * h =  (TH2D*) h3D->Project3D("zxe");
  if(scaleWidth) h ->Scale(1.,"width");

  return h;

}


TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id, 
						       Float_t minPt, Float_t maxPt,
						       Float_t minEta, Float_t maxEta,
						       Bool_t scaleWidth) { 
  // Returns a projection of the 3D histo pt/eta/vz.
  // WARNING: since that is a histo, the requested range will be discretized to the binning.
  // Always avoids under (over) flows
  // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)

  TH3D * h3D = GetHistoPtEtaVz(id);

  // Get range in terms of bin numners.  If the float range is
  // less than -11111 take the range from the first to the last bin (i.e. no
  // under/over-flows)
  Int_t min1  = minPt  < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
  Int_t min2  = minEta < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minEta);

  Int_t max1  = maxPt  < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
  Int_t max2  = maxEta < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);


  TString hname = h3D->GetName();
  hname = hname +  "_Vz_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);

  if (gROOT->FindObjectAny(hname.Data())){
    AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
    hname+="_2";
  }

  TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E");
  if(scaleWidth) h ->Scale(1.,"width");
  return h;


}

TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id, 
							Float_t minPt, Float_t maxPt, 
							Float_t minVz, Float_t maxVz,
							Bool_t scaleWidth) {
  // Returns a projection of the 3D histo pt/eta/vz.
  // WARNING: since that is a histo, the requested range will be discretized to the binning.
  // Always avoids under (over) flows
  // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)

  TH3D * h3D = GetHistoPtEtaVz(id);

  // Get range in terms of bin numners.  If the float range is
  // less than -11111 take the range from the first to the last bin (i.e. no
  // under/over-flows)
  Int_t min1 = minPt < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
  Int_t min2 = minVz < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minVz);

  Int_t max1 = maxPt < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
  Int_t max2 = maxVz < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxVz-0.00001);

  TString hname = h3D->GetName();
  hname = hname +  "_Eta_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);

  if (gROOT->FindObjectAny(hname.Data())){
    AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
    hname+="_2";
  }

  TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E");
  if(scaleWidth) h ->Scale(1.,"width");
  return h;
}

TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoProcess(Histo_t id) {

  // Returns histogram with particle specties

  TString name = TString(kHistoPrefix[id])+"_Process";

  TH1D * h =  (TH1D*) GetHisto(name);
  if (!h) {
    name+=fHNameSuffix;
    Bool_t oldStatus = TH1::AddDirectoryStatus();
    TH1::AddDirectory(kFALSE);

    AliInfo(Form("Booking histo %s",name.Data()));

    h = new TH1D (name.Data(), Form("Particle production process (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);			 
    Int_t nbin = kPNoProcess+1;
    for(Int_t ibin = 0; ibin < nbin; ibin++){
      h->GetXaxis()->SetBinLabel(ibin+1,TMCProcessName[ibin]);      
    }
    TH1::AddDirectory(oldStatus);
    fList->Add(h);


  }
  return h;
  

}


TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoSpecies(Histo_t id) {

  // Returns histogram with particle specties

  TString name = TString(kHistoPrefix[id])+"_Species";

  TH1D * h =  (TH1D*) GetHisto(name);
  if (!h) {
    name+=fHNameSuffix;
    Bool_t oldStatus = TH1::AddDirectoryStatus();
    TH1::AddDirectory(kFALSE);

    AliInfo(Form("Booking histo %s",name.Data()));

    h = new TH1D (name.Data(), Form("Particle species (%s)",kHistoPrefix[id]), kNPart+1, -0.5, kNPart+1-0.5);			 
    Int_t nbin = kNPart+1;
    for(Int_t ibin = 0; ibin < nbin; ibin++){
      h->GetXaxis()->SetBinLabel(ibin+1,kSpeciesName[ibin]);      
    }
    TH1::AddDirectory(oldStatus);
    fList->Add(h);


  }
  return h;
  

}

TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVzEvent(Histo_t id) {

  // Returns histogram with Vz of the event

  TString name = TString(kHistoPrefix[id])+"_VzEvent";

  TH1D * h =  (TH1D*) GetHisto(name);
  if (!h) {
    name+=fHNameSuffix;
    Bool_t oldStatus = TH1::AddDirectoryStatus();
    TH1::AddDirectory(kFALSE);

    AliInfo(Form("Booking histo %s",name.Data()));

    h = new TH1D (name.Data(), Form("Vz of the event (%s)",kHistoPrefix[id]), 10, -10, 10);			 
    h->Sumw2();
    h->SetXTitle("V_{z}^{event} (cm)");
    TH1::AddDirectory(oldStatus);
    fList->Add(h);


  }
  return h;
  

}


TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() {
  // Returns histogram with event statistiscs (processed events at each step)

  TH1I * h =  (TH1I*) GetHisto("hStats");
  if (!h) h = BookHistoStats();
  return h;
  
} 


TH1 * AliAnalysisMultPbTrackHistoManager::GetHisto(const char * name) {
  //Search list for histo
  // TODO: keep track of histo id rather than searching by name?
  return (TH1*) fList->FindObject(TString(name)+fHNameSuffix);

}

void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * option) {
  // Scales all histos in the list for nev
  // option can be used to pass further options to TH1::Scale
  TH1 * h = 0;
  TIter iter = fList->MakeIterator();
  while ((h = (TH1*) iter.Next())) {
    if (!h->InheritsFrom("TH1")) {
      AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName()));
    }
    AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
    
    h->Scale(1./nev,option);
  }

}

TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, const char * title) {
  // Books a 3D histo of Pt/eta/vtx
  // TODO: make the binning settable, variable binning?

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  TString hname = name;
  hname+=fHNameSuffix;

  AliInfo(Form("Booking %s",hname.Data()));
  
  // Binning from Jacek task
  // const Int_t nptbins = 49;
  // const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0};//,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
  const Int_t nptbins = 68;
  const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
  
  const Int_t netabins=20;
  Double_t binsEta[netabins+1];
  Float_t minEta = -1;
  Float_t maxEta =  1;
  Float_t etaStep = (maxEta-minEta)/netabins;
  for(Int_t ibin = 0; ibin < netabins; ibin++){    
    binsEta[ibin]   = minEta + ibin*etaStep;
    binsEta[ibin+1] = minEta + ibin*etaStep + etaStep;
  }

  const Int_t nvzbins=10;
  Double_t binsVz[nvzbins+1];
  Float_t minVz = -10;
  Float_t maxVz =  10;
  Float_t vzStep = (maxVz-minVz)/nvzbins;
  for(Int_t ibin = 0; ibin < nvzbins; ibin++){    
    binsVz[ibin]   = minVz + ibin*vzStep;
    binsVz[ibin+1] = minVz + ibin*vzStep + vzStep;
  }
 

  TH3D * h = new TH3D (hname,title, 
		       nptbins,  binsPt,
		       netabins, binsEta,
		       nvzbins,  binsVz
		       );

  h->SetXTitle("p_{T} (GeV)");
  h->SetYTitle("#eta");
  h->SetZTitle("V_{z}^{tracks} (cm)");
  h->Sumw2();
  
  fList->Add(h);

  TH1::AddDirectory(oldStatus);
  return h;
}

TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
  // Books a DCA histo 

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  TString hname = name;
  hname+=fHNameSuffix;

  AliInfo(Form("Booking %s",hname.Data()));
  

  TH1D * h = new TH1D (hname,title, 200,0,200);

  h->SetXTitle("#Delta DCA");
  h->Sumw2();
  
  fList->Add(h);

  TH1::AddDirectory(oldStatus);
  return h;
}
TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoMult(const char * name, const char * title) {
  // Books a multiplicity histo 

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  TString hname = name;
  hname+=fHNameSuffix;

  AliInfo(Form("Booking %s",hname.Data()));
  

  TH1D * h = new TH1D (hname,title, 600, 0,6000);

  h->SetXTitle("N tracks");
  h->Sumw2();
  
  fList->Add(h);

  TH1::AddDirectory(oldStatus);
  return h;
}

TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() {
  // Books histogram with event statistiscs (processed events at each step)

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  AliInfo(Form("Booking Stat histo"));

  TH1I * h = new TH1I (TString("hStats")+fHNameSuffix, "Number of processed events", kNStatBins, -0.5, kNStatBins-0.5);
  for(Int_t istatbin = 0; istatbin < kNStatBins; istatbin++){
    h->GetXaxis()->SetBinLabel(istatbin+1,kStatStepNames[istatbin]);
  }
  TH1::AddDirectory(oldStatus);
  fList->Add(h);
  return h;
}

Int_t AliAnalysisMultPbTrackHistoManager::GetLocalParticleID(AliMCParticle * part) {
  // returns the local code (Part_t)
  
  Int_t pdgcode = part->PdgCode();
  switch(pdgcode) {

  case 211:
    return kPartPiPlus;
    break;
  case -211:
    return kPartPiMinus;
    break;
  case 2212:
    return kPartP;
    break;
  case -2212:
    return kPartPBar;
    break;
  case 321:
    return kPartKPlus;
    break;
  case -321:
    return kPartKMinus;
    break;
  case -11:
    return kPartLMinus;
    break;
  case 11:
    return kPartLPlus;
    break;
  case -13:
    return kPartLMinus;
    break;
  case 13:
    return kPartLPlus;
    break;
  default:
    return kPartOther;
  }
}
Commit	Line	Data
a23f7c97	1	#include "AliAnalysisMultPbTrackHistoManager.h"
	2	#include "AliLog.h"
	3	#include "TH1.h"
	4	#include "TH3D.h"
	5	#include "TH1I.h"
	6	#include "TROOT.h"
e0376287	7	#include "TMCProcess.h"
3b8cbf2d	8	#include "AliMCParticle.h"
e0376287	9
a23f7c97	10	#include <iostream>
aa6151eb	11	#include "TH2D.h"
3b8cbf2d	12
a23f7c97	13	using namespace std;
3b8cbf2d	14
a23f7c97	15	ClassImp(AliAnalysisMultPbTrackHistoManager)
a23f7c97	16
8b628a20	17	const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[] = { "All Events", "After physics Selection", "After centrality selection", "With Vertex" };
a23f7c97	18	const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim",
	19	"hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
	20	const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[] = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
e0376287	21	const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[] = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake"};
4d0aa70f	22	const char * AliAnalysisMultPbTrackHistoManager::kSpeciesName[] = { "pi+", "K+", "p", "l+", "pi-", "K-", "barp", "l-", "Others"};
a23f7c97	23
	24
	25	AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager() : AliHistoListWrapper(), fHNameSuffix(""){
	26	// standard ctor
	27
	28	}
	29
	30	AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const char * name, const char * title): AliHistoListWrapper(name,title), fHNameSuffix("") {
	31	//named ctor
	32	};
	33
	34	AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) : AliHistoListWrapper (obj) {
	35	// copy ctor
	36	AliError("Not Implemented");
	37	};
	38
	39	AliAnalysisMultPbTrackHistoManager::~AliAnalysisMultPbTrackHistoManager() {
	40	// dtor
	41
	42	}
	43
4d0aa70f	44	TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id, Int_t particle) {
a23f7c97	45	// Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
a23f7c97	46
4d0aa70f	47	TString name = kHistoPtEtaVzNames[id];
	48
	49	if(particle >= 0) name += kSpeciesName[particle];
	50
	51	TH3D * h = (TH3D*) GetHisto(name.Data());
a23f7c97	52	if (!h) {
4d0aa70f	53	h = BookHistoPtEtaVz(name.Data(), Form("Pt Eta Vz distribution (%s)",kHistoPtEtaVzNames[id]));
a23f7c97	54	}
	55
	56	return h;
	57
	58	}
	59
	60	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
	61	// Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
	62
	63	TH1D * h = (TH1D*) GetHisto(kHistoDCANames[id]);
	64	if (!h) {
	65	h = BookHistoDCA(kHistoDCANames[id], Form("Pt Eta Vz distribution (%s)",kHistoDCANames[id]));
	66	}
	67
	68	return h;
	69
	70	}
	71
e0376287	72	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMult(Histo_t id) {
	73	// Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
	74
	75	TString name = kHistoPrefix[id];
	76	name += "Mult";
	77	TH1D * h = (TH1D*) GetHisto(name.Data());
	78	if (!h) {
	79	h = BookHistoMult(name.Data(), Form("Multiplicity distribution (%s)",kHistoPrefix[id]));
	80	}
	81
	82	return h;
	83
	84	}
	85
a23f7c97	86	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id,
	87	Float_t minEta, Float_t maxEta,
	88	Float_t minVz, Float_t maxVz,
	89	Bool_t scaleWidth) {
	90	// Returns a projection of the 3D histo pt/eta/vz.
	91	// WARNING: since that is a histo, the requested range will be discretized to the binning.
	92	// Always avoids under (over) flows
	93	// If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
	94
	95	TH3D * h3D = GetHistoPtEtaVz(id);
	96
	97	// Get range in terms of bin numners. If the float range is
	98	// less than -11111 take the range from the first to the last bin (i.e. no
	99	// under/over-flows)
4d0aa70f	100	// FIXME: UNDERFLOWS
	101	Int_t min1 = minEta < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
	102	Int_t min2 = minVz < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;
a23f7c97	103
4d0aa70f	104	// Int_t max1 = maxEta < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
	105	// Int_t max2 = maxVz < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
	106	Int_t max1 = maxEta < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
	107	Int_t max2 = maxVz < -11111 ? -1 : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
a23f7c97	108
	109
	110	TString hname = h3D->GetName();
	111	hname = hname + "_pt_" + long (min1) + "_" + long(max1) + "_" + long (min2) + "_" + long(max2);
	112
	113
	114	if (gROOT->FindObjectAny(hname.Data())){
e0376287	115	AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
e0376287	116	hname += "_2";
a23f7c97	117	}
	118
	119	TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E");
	120	if(scaleWidth) h ->Scale(1.,"width");
	121
	122	return h;
	123
	124	}
	125
aa6151eb	126
	127	TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoPtVz (Histo_t id,
	128	Float_t minEta, Float_t maxEta,
	129	Bool_t scaleWidth) {
	130	// Returns a projection of the 3D histo pt/eta/vz.
	131	// WARNING: since that is a histo, the requested range will be discretized to the binning.
	132	// Always avoids under (over) flows
	133	// If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
	134
	135	// FIXME: what do I do here for the scaling?
	136
	137	TH3D * h3D = GetHistoPtEtaVz(id);
	138
	139	// Get range in terms of bin numners. If the float range is
	140	// less than -11111 take the range from the first to the last bin (i.e. no
	141	// under/over-flows)
	142	Int_t min1 = minEta < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
	143	Int_t max1 = maxEta < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
	144
	145
	146	TString hname = h3D->GetName();
	147	hname = hname + "_ptvz_" + long (min1) + "_" + long(max1);
	148
	149	if (gROOT->FindObjectAny(hname.Data())){
	150	AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
	151	hname += "_2";
	152	}
	153
	154	h3D->GetYaxis()->SetRange(min1,max1);
	155
	156	TH2D * h = (TH2D*) h3D->Project3D("zxe");
	157	if(scaleWidth) h ->Scale(1.,"width");
	158
	159	return h;
	160
	161	}
	162
	163
a23f7c97	164	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id,
	165	Float_t minPt, Float_t maxPt,
	166	Float_t minEta, Float_t maxEta,
	167	Bool_t scaleWidth) {
	168	// Returns a projection of the 3D histo pt/eta/vz.
	169	// WARNING: since that is a histo, the requested range will be discretized to the binning.
	170	// Always avoids under (over) flows
	171	// If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
	172
	173	TH3D * h3D = GetHistoPtEtaVz(id);
	174
	175	// Get range in terms of bin numners. If the float range is
	176	// less than -11111 take the range from the first to the last bin (i.e. no
	177	// under/over-flows)
	178	Int_t min1 = minPt < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
	179	Int_t min2 = minEta < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minEta);
	180
	181	Int_t max1 = maxPt < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
	182	Int_t max2 = maxEta < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
	183
	184
	185	TString hname = h3D->GetName();
	186	hname = hname + "_Vz_" + long (min1) + "_" + long(max1) + "_" + long (min2) + "_" + long(max2);
	187
	188	if (gROOT->FindObjectAny(hname.Data())){
e0376287	189	AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
e0376287	190	hname+="_2";
a23f7c97	191	}
	192
	193	TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E");
	194	if(scaleWidth) h ->Scale(1.,"width");
	195	return h;
	196
	197
	198	}
	199
	200	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id,
	201	Float_t minPt, Float_t maxPt,
	202	Float_t minVz, Float_t maxVz,
	203	Bool_t scaleWidth) {
	204	// Returns a projection of the 3D histo pt/eta/vz.
	205	// WARNING: since that is a histo, the requested range will be discretized to the binning.
	206	// Always avoids under (over) flows
	207	// If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
	208
	209	TH3D * h3D = GetHistoPtEtaVz(id);
	210
	211	// Get range in terms of bin numners. If the float range is
	212	// less than -11111 take the range from the first to the last bin (i.e. no
	213	// under/over-flows)
	214	Int_t min1 = minPt < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
	215	Int_t min2 = minVz < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minVz);
	216
	217	Int_t max1 = maxPt < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
	218	Int_t max2 = maxVz < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxVz-0.00001);
	219
	220	TString hname = h3D->GetName();
	221	hname = hname + "_Eta_" + long (min1) + "_" + long(max1) + "_" + long (min2) + "_" + long(max2);
	222
	223	if (gROOT->FindObjectAny(hname.Data())){
e0376287	224	AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
e0376287	225	hname+="_2";
a23f7c97	226	}
	227
	228	TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E");
	229	if(scaleWidth) h ->Scale(1.,"width");
	230	return h;
	231	}
	232
3b8cbf2d	233	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoProcess(Histo_t id) {
e0376287	234
	235	// Returns histogram with particle specties
	236
3b8cbf2d	237	TString name = TString(kHistoPrefix[id])+"_Process";
e0376287	238
	239	TH1D * h = (TH1D*) GetHisto(name);
	240	if (!h) {
	241	name+=fHNameSuffix;
	242	Bool_t oldStatus = TH1::AddDirectoryStatus();
	243	TH1::AddDirectory(kFALSE);
	244
	245	AliInfo(Form("Booking histo %s",name.Data()));
	246
3b8cbf2d	247	h = new TH1D (name.Data(), Form("Particle production process (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);
e0376287	248	Int_t nbin = kPNoProcess+1;
	249	for(Int_t ibin = 0; ibin < nbin; ibin++){
	250	h->GetXaxis()->SetBinLabel(ibin+1,TMCProcessName[ibin]);
	251	}
	252	TH1::AddDirectory(oldStatus);
	253	fList->Add(h);
	254
	255
	256	}
	257	return h;
	258
	259
	260	}
	261
	262
3b8cbf2d	263	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoSpecies(Histo_t id) {
	264
	265	// Returns histogram with particle specties
	266
	267	TString name = TString(kHistoPrefix[id])+"_Species";
	268
	269	TH1D * h = (TH1D*) GetHisto(name);
	270	if (!h) {
	271	name+=fHNameSuffix;
	272	Bool_t oldStatus = TH1::AddDirectoryStatus();
	273	TH1::AddDirectory(kFALSE);
	274
	275	AliInfo(Form("Booking histo %s",name.Data()));
	276
	277	h = new TH1D (name.Data(), Form("Particle species (%s)",kHistoPrefix[id]), kNPart+1, -0.5, kNPart+1-0.5);
	278	Int_t nbin = kNPart+1;
	279	for(Int_t ibin = 0; ibin < nbin; ibin++){
	280	h->GetXaxis()->SetBinLabel(ibin+1,kSpeciesName[ibin]);
	281	}
	282	TH1::AddDirectory(oldStatus);
	283	fList->Add(h);
	284
	285
	286	}
	287	return h;
	288
	289
	290	}
	291
4d0aa70f	292	TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVzEvent(Histo_t id) {
	293
	294	// Returns histogram with Vz of the event
	295
	296	TString name = TString(kHistoPrefix[id])+"_VzEvent";
	297
	298	TH1D * h = (TH1D*) GetHisto(name);
	299	if (!h) {
	300	name+=fHNameSuffix;
	301	Bool_t oldStatus = TH1::AddDirectoryStatus();
	302	TH1::AddDirectory(kFALSE);
	303
	304	AliInfo(Form("Booking histo %s",name.Data()));
	305
	306	h = new TH1D (name.Data(), Form("Vz of the event (%s)",kHistoPrefix[id]), 10, -10, 10);
	307	h->Sumw2();
	308	h->SetXTitle("V_{z}^{event} (cm)");
	309	TH1::AddDirectory(oldStatus);
	310	fList->Add(h);
	311
	312
	313	}
	314	return h;
	315
	316
	317	}
	318
3b8cbf2d	319
a23f7c97	320
	321	TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() {
	322	// Returns histogram with event statistiscs (processed events at each step)
	323
	324	TH1I * h = (TH1I*) GetHisto("hStats");
	325	if (!h) h = BookHistoStats();
	326	return h;
	327
	328	}
	329
	330
	331
	332	TH1 * AliAnalysisMultPbTrackHistoManager::GetHisto(const char * name) {
	333	//Search list for histo
	334	// TODO: keep track of histo id rather than searching by name?
	335	return (TH1*) fList->FindObject(TString(name)+fHNameSuffix);
	336
	337	}
	338
	339	void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * option) {
	340	// Scales all histos in the list for nev
	341	// option can be used to pass further options to TH1::Scale
	342	TH1 * h = 0;
	343	TIter iter = fList->MakeIterator();
	344	while ((h = (TH1*) iter.Next())) {
	345	if (!h->InheritsFrom("TH1")) {
	346	AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName()));
	347	}
e0376287	348	AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
a23f7c97	349
	350	h->Scale(1./nev,option);
	351	}
	352
	353	}
	354
	355	TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, const char * title) {
	356	// Books a 3D histo of Pt/eta/vtx
	357	// TODO: make the binning settable, variable binning?
	358
	359	Bool_t oldStatus = TH1::AddDirectoryStatus();
	360	TH1::AddDirectory(kFALSE);
	361
	362	TString hname = name;
	363	hname+=fHNameSuffix;
	364
	365	AliInfo(Form("Booking %s",hname.Data()));
	366
e0376287	367	// Binning from Jacek task
eef42d18	368	// const Int_t nptbins = 49;
eef42d18	369	// const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0};//,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
e0376287	370	const Int_t nptbins = 68;
	371	const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
	372
	373	const Int_t netabins=20;
	374	Double_t binsEta[netabins+1];
	375	Float_t minEta = -1;
	376	Float_t maxEta = 1;
	377	Float_t etaStep = (maxEta-minEta)/netabins;
	378	for(Int_t ibin = 0; ibin < netabins; ibin++){
	379	binsEta[ibin] = minEta + ibin*etaStep;
	380	binsEta[ibin+1] = minEta + ibin*etaStep + etaStep;
	381	}
	382
	383	const Int_t nvzbins=10;
	384	Double_t binsVz[nvzbins+1];
	385	Float_t minVz = -10;
	386	Float_t maxVz = 10;
	387	Float_t vzStep = (maxVz-minVz)/nvzbins;
	388	for(Int_t ibin = 0; ibin < nvzbins; ibin++){
	389	binsVz[ibin] = minVz + ibin*vzStep;
	390	binsVz[ibin+1] = minVz + ibin*vzStep + vzStep;
	391	}
	392
a23f7c97	393
a23f7c97	394	TH3D * h = new TH3D (hname,title,
e0376287	395	nptbins, binsPt,
	396	netabins, binsEta,
	397	nvzbins, binsVz
a23f7c97	398	);
a23f7c97	399
4d0aa70f	400	h->SetXTitle("p_{T} (GeV)");
e0376287	401	h->SetYTitle("#eta");
4d0aa70f	402	h->SetZTitle("V_{z}^{tracks} (cm)");
a23f7c97	403	h->Sumw2();
	404
	405	fList->Add(h);
	406
	407	TH1::AddDirectory(oldStatus);
	408	return h;
	409	}
	410
	411	TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
	412	// Books a DCA histo
	413
	414	Bool_t oldStatus = TH1::AddDirectoryStatus();
	415	TH1::AddDirectory(kFALSE);
	416
	417	TString hname = name;
	418	hname+=fHNameSuffix;
	419
	420	AliInfo(Form("Booking %s",hname.Data()));
	421
	422
54b31d6a	423	TH1D * h = new TH1D (hname,title, 200,0,200);
a23f7c97	424
	425	h->SetXTitle("#Delta DCA");
	426	h->Sumw2();
	427
	428	fList->Add(h);
	429
	430	TH1::AddDirectory(oldStatus);
	431	return h;
	432	}
e0376287	433	TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoMult(const char * name, const char * title) {
	434	// Books a multiplicity histo
	435
	436	Bool_t oldStatus = TH1::AddDirectoryStatus();
	437	TH1::AddDirectory(kFALSE);
	438
	439	TString hname = name;
	440	hname+=fHNameSuffix;
	441
	442	AliInfo(Form("Booking %s",hname.Data()));
	443
	444
	445	TH1D * h = new TH1D (hname,title, 600, 0,6000);
	446
	447	h->SetXTitle("N tracks");
	448	h->Sumw2();
	449
	450	fList->Add(h);
	451
	452	TH1::AddDirectory(oldStatus);
	453	return h;
	454	}
a23f7c97	455
	456	TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() {
	457	// Books histogram with event statistiscs (processed events at each step)
	458
	459	Bool_t oldStatus = TH1::AddDirectoryStatus();
	460	TH1::AddDirectory(kFALSE);
	461
	462	AliInfo(Form("Booking Stat histo"));
	463
	464	TH1I * h = new TH1I (TString("hStats")+fHNameSuffix, "Number of processed events", kNStatBins, -0.5, kNStatBins-0.5);
	465	for(Int_t istatbin = 0; istatbin < kNStatBins; istatbin++){
	466	h->GetXaxis()->SetBinLabel(istatbin+1,kStatStepNames[istatbin]);
	467	}
	468	TH1::AddDirectory(oldStatus);
	469	fList->Add(h);
	470	return h;
	471	}
	472
3b8cbf2d	473	Int_t AliAnalysisMultPbTrackHistoManager::GetLocalParticleID(AliMCParticle * part) {
	474	// returns the local code (Part_t)
	475
	476	Int_t pdgcode = part->PdgCode();
	477	switch(pdgcode) {
	478
	479	case 211:
	480	return kPartPiPlus;
	481	break;
	482	case -211:
	483	return kPartPiMinus;
	484	break;
	485	case 2212:
	486	return kPartP;
	487	break;
	488	case -2212:
	489	return kPartPBar;
	490	break;
	491	case 321:
	492	return kPartKPlus;
	493	break;
	494	case -321:
	495	return kPartKMinus;
	496	break;
	497	case -11:
	498	return kPartLMinus;
	499	break;
	500	case 11:
	501	return kPartLPlus;
	502	break;
	503	case -13:
	504	return kPartLMinus;
	505	break;
	506	case 13:
	507	return kPartLPlus;
	508	break;
	509	default:
	510	return kPartOther;
	511	}
	512	}
a23f7c97	513
	514
	515