[u/mrichter/AliRoot.git] / PWG0 / AliPWG0Helper.cxx

/* $Id$ */

#include <AliPWG0Helper.h>

#include <TParticle.h>
#include <TParticlePDG.h>
#include <TH1.h>
#include <TH3.h>

#include <AliLog.h>
#include <AliESD.h>
#include <AliESDVertex.h>

//____________________________________________________________________
ClassImp(AliPWG0Helper)

//____________________________________________________________________
Bool_t AliPWG0Helper::IsEventTriggered(AliESD* aEsd)
{
  // check if the event was triggered
  //
  // this function needs the branch fTriggerMask
  //
  // MB should be
  // ITS_SPD_GFO_L0  : 32
  // VZERO_OR_LEFT   : 1
  // VZERO_OR_RIGHT  : 2

  ULong64_t triggerMask = aEsd->GetTriggerMask();

  if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2)))
    return kTRUE;

  return kFALSE;
}

//____________________________________________________________________
Bool_t AliPWG0Helper::IsVertexReconstructed(AliESD* aEsd)
{
  // checks if the vertex is reasonable
  //
  // this function needs the branches fSPDVertex*


  const AliESDVertex* vtxESD = aEsd->GetVertex();

  // the vertex should be reconstructed
  if (strcmp(vtxESD->GetName(), "default")==0)
    return kFALSE;

  Double_t vtx_res[3];
  vtx_res[0] = vtxESD->GetXRes();
  vtx_res[1] = vtxESD->GetYRes();
  vtx_res[2] = vtxESD->GetZRes();

  if (vtx_res[2]==0 || vtx_res[2]>0.1)
    return kFALSE;

  // check Ncontributors, if <0 it means error *gna*

  return kTRUE;
}

//____________________________________________________________________
Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
{
  //
  // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
  // shall be counted as a primary particle
  //
  // This function or a equivalent should be available in some common place of AliRoot
  //

  // if the particle has a daughter primary, we do not want to count it
  if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
  {
    if (adebug)
      printf("Dropping particle because it has a daughter among the primaries.\n");
    return kFALSE;
  }

  Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());

  // skip quarks and gluon
  if (pdgCode <= 10 || pdgCode == 21)
  {
    if (adebug)
      printf("Dropping particle because it is a quark or gluon.\n");
    return kFALSE;
  }

  if (strcmp(aParticle->GetName(),"XXX") == 0)
  {
    if (adebug)
      printf("WARNING: There is a particle named XXX.\n");
    return kFALSE;
  }

  TParticlePDG* pdgPart = aParticle->GetPDG();

  if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
  {
    if (adebug)
      printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode);
    return kFALSE;
  }

  if (pdgPart->Charge() == 0)
  {
    if (adebug)
      printf("Dropping particle because it is not charged.\n");
    return kFALSE;
  }

  return kTRUE;
}

//____________________________________________________________________
void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save)
{
  // create projections of 3d hists to all 2d combinations
  // the histograms are not returned, just use them from memory or use this to create them in a file

  TH1* proj = hist->Project3D("yx");
  proj->SetXTitle(hist->GetXaxis()->GetTitle());
  proj->SetYTitle(hist->GetYaxis()->GetTitle());
  if (save)
    proj->Write();

  proj = hist->Project3D("zx");
  proj->SetXTitle(hist->GetXaxis()->GetTitle());
  proj->SetYTitle(hist->GetZaxis()->GetTitle());
  if (save)
    proj->Write();

  proj = hist->Project3D("zy");
  proj->SetXTitle(hist->GetYaxis()->GetTitle());
  proj->SetYTitle(hist->GetZaxis()->GetTitle());
  if (save)
    proj->Write();
}

//____________________________________________________________________
void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save)
{
  // create projections of the 3d hists divides them
  // axis decides to which plane, if axis is 0 to all planes
  // the histograms are not returned, just use them from memory or use this to create them in a file

  if (axis == 0)
  {
    CreateDividedProjections(hist, hist2, "yx", putErrors, save);
    CreateDividedProjections(hist, hist2, "zx", putErrors, save);
    CreateDividedProjections(hist, hist2, "zy", putErrors, save);

    return;
  }

  TH1* proj = hist->Project3D(axis);

  if (strlen(axis) == 2)
  {
    proj->SetYTitle(GetAxisTitle(hist, axis[0]));
    proj->SetXTitle(GetAxisTitle(hist, axis[1]));
  }
  else if (strlen(axis) == 1)
    proj->SetXTitle(GetAxisTitle(hist, axis[0]));

  TH1* proj2 = hist2->Project3D(axis);
  if (strlen(axis) == 2)
  {
    proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
    proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
  }
  else if (strlen(axis) == 1)
    proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));

  TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
  //printf("doing axis: %s, x axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsX(), proj2->GetNbinsX(), division->GetXaxis()->GetBinLowEdge(1), proj2->GetXaxis()->GetBinLowEdge(1), division->GetXaxis()->GetBinUpEdge(division->GetNbinsX()), proj2->GetXaxis()->GetBinUpEdge(proj2->GetNbinsX()));
  //printf("doing axis: %s, y axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsY(), proj2->GetNbinsY(), division->GetYaxis()->GetBinLowEdge(1), proj2->GetYaxis()->GetBinLowEdge(1), division->GetYaxis()->GetBinUpEdge(division->GetNbinsY()), proj2->GetYaxis()->GetBinUpEdge(proj2->GetNbinsY()));
  division->Divide(proj, proj2, 1, 1, "B");
  division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle()));

  if (putErrors)
  {
    division->Sumw2();
    if (division->GetDimension() == 1)
    {
      Int_t nBins = division->GetNbinsX();
      for (Int_t i = 1; i <= nBins; ++i)
        if (proj2->GetBinContent(i) != 0)
          division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i));
    }
    else if (division->GetDimension() == 2)
    {
      Int_t nBinsX = division->GetNbinsX();
      Int_t nBinsY = division->GetNbinsY();
      for (Int_t i = 1; i <= nBinsX; ++i)
        for (Int_t j = 1; j <= nBinsY; ++j)
          if (proj2->GetBinContent(i, j) != 0)
            division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
    }
  }

  if (save)
  {
    proj->Write();
    proj2->Write();
    division->Write();
  }
}

//____________________________________________________________________
const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
{
  // returns the title of the axis given in axis (x, y, z)

  if (axis == 'x')
    return hist->GetXaxis()->GetTitle();
  else if (axis == 'y')
    return hist->GetYaxis()->GetTitle();
  else if (axis == 'z')
    return hist->GetZaxis()->GetTitle();

  return 0;
}
Commit	Line	Data
04a7657f	1	/* $Id$ */
	2
	3	#include <AliPWG0Helper.h>
	4
	5	#include <TParticle.h>
	6	#include <TParticlePDG.h>
847489f7	7	#include <TH1.h>
25db2d85	8	#include <TH3.h>
04a7657f	9
	10	#include <AliLog.h>
	11	#include <AliESD.h>
	12	#include <AliESDVertex.h>
	13
	14	//____________________________________________________________________
	15	ClassImp(AliPWG0Helper)
	16
	17	//____________________________________________________________________
	18	Bool_t AliPWG0Helper::IsEventTriggered(AliESD* aEsd)
	19	{
	20	// check if the event was triggered
	21	//
5c495d37	22	// this function needs the branch fTriggerMask
5c495d37	23	//
04a7657f	24	// MB should be
	25	// ITS_SPD_GFO_L0 : 32
	26	// VZERO_OR_LEFT : 1
	27	// VZERO_OR_RIGHT : 2
	28
	29	ULong64_t triggerMask = aEsd->GetTriggerMask();
	30
	31	if (triggerMask&32 && ((triggerMask&1) \|\| (triggerMask&2)))
	32	return kTRUE;
	33
	34	return kFALSE;
	35	}
	36
	37	//____________________________________________________________________
	38	Bool_t AliPWG0Helper::IsVertexReconstructed(AliESD* aEsd)
	39	{
	40	// checks if the vertex is reasonable
5c495d37	41	//
	42	// this function needs the branches fSPDVertex*
	43
04a7657f	44
	45	const AliESDVertex* vtxESD = aEsd->GetVertex();
	46
	47	// the vertex should be reconstructed
	48	if (strcmp(vtxESD->GetName(), "default")==0)
	49	return kFALSE;
	50
	51	Double_t vtx_res[3];
	52	vtx_res[0] = vtxESD->GetXRes();
	53	vtx_res[1] = vtxESD->GetYRes();
	54	vtx_res[2] = vtxESD->GetZRes();
	55
	56	if (vtx_res[2]==0 \|\| vtx_res[2]>0.1)
	57	return kFALSE;
	58
9e952c39	59	// check Ncontributors, if <0 it means error gna
9e952c39	60
04a7657f	61	return kTRUE;
	62	}
	63
	64	//____________________________________________________________________
7584d357	65	Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
04a7657f	66	{
04a7657f	67	//
25db2d85	68	// this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
	69	// shall be counted as a primary particle
	70	//
04a7657f	71	// This function or a equivalent should be available in some common place of AliRoot
	72	//
	73
	74	// if the particle has a daughter primary, we do not want to count it
	75	if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
	76	{
7584d357	77	if (adebug)
25db2d85	78	printf("Dropping particle because it has a daughter among the primaries.\n");
04a7657f	79	return kFALSE;
	80	}
	81
	82	Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
	83
	84	// skip quarks and gluon
	85	if (pdgCode <= 10 \|\| pdgCode == 21)
	86	{
7584d357	87	if (adebug)
25db2d85	88	printf("Dropping particle because it is a quark or gluon.\n");
04a7657f	89	return kFALSE;
	90	}
	91
	92	if (strcmp(aParticle->GetName(),"XXX") == 0)
	93	{
7584d357	94	if (adebug)
25db2d85	95	printf("WARNING: There is a particle named XXX.\n");
04a7657f	96	return kFALSE;
	97	}
	98
	99	TParticlePDG* pdgPart = aParticle->GetPDG();
	100
	101	if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
	102	{
7584d357	103	if (adebug)
25db2d85	104	printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode);
04a7657f	105	return kFALSE;
	106	}
	107
	108	if (pdgPart->Charge() == 0)
	109	{
7584d357	110	if (adebug)
25db2d85	111	printf("Dropping particle because it is not charged.\n");
04a7657f	112	return kFALSE;
	113	}
	114
	115	return kTRUE;
	116	}
847489f7	117
847489f7	118	//____________________________________________________________________
29771dc8	119	void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save)
847489f7	120	{
	121	// create projections of 3d hists to all 2d combinations
	122	// the histograms are not returned, just use them from memory or use this to create them in a file
	123
	124	TH1* proj = hist->Project3D("yx");
	125	proj->SetXTitle(hist->GetXaxis()->GetTitle());
	126	proj->SetYTitle(hist->GetYaxis()->GetTitle());
29771dc8	127	if (save)
29771dc8	128	proj->Write();
847489f7	129
	130	proj = hist->Project3D("zx");
	131	proj->SetXTitle(hist->GetXaxis()->GetTitle());
	132	proj->SetYTitle(hist->GetZaxis()->GetTitle());
29771dc8	133	if (save)
29771dc8	134	proj->Write();
847489f7	135
	136	proj = hist->Project3D("zy");
	137	proj->SetXTitle(hist->GetYaxis()->GetTitle());
	138	proj->SetYTitle(hist->GetZaxis()->GetTitle());
29771dc8	139	if (save)
29771dc8	140	proj->Write();
847489f7	141	}
1afae8ff	142
1afae8ff	143	//____________________________________________________________________
29771dc8	144	void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save)
1afae8ff	145	{
	146	// create projections of the 3d hists divides them
	147	// axis decides to which plane, if axis is 0 to all planes
	148	// the histograms are not returned, just use them from memory or use this to create them in a file
	149
	150	if (axis == 0)
	151	{
29771dc8	152	CreateDividedProjections(hist, hist2, "yx", putErrors, save);
	153	CreateDividedProjections(hist, hist2, "zx", putErrors, save);
	154	CreateDividedProjections(hist, hist2, "zy", putErrors, save);
1afae8ff	155
	156	return;
	157	}
	158
	159	TH1* proj = hist->Project3D(axis);
0ab29cfa	160
	161	if (strlen(axis) == 2)
	162	{
	163	proj->SetYTitle(GetAxisTitle(hist, axis[0]));
	164	proj->SetXTitle(GetAxisTitle(hist, axis[1]));
	165	}
	166	else if (strlen(axis) == 1)
	167	proj->SetXTitle(GetAxisTitle(hist, axis[0]));
1afae8ff	168
1afae8ff	169	TH1* proj2 = hist2->Project3D(axis);
0ab29cfa	170	if (strlen(axis) == 2)
	171	{
	172	proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
	173	proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
	174	}
	175	else if (strlen(axis) == 1)
	176	proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));
1afae8ff	177
1afae8ff	178	TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
29771dc8	179	//printf("doing axis: %s, x axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsX(), proj2->GetNbinsX(), division->GetXaxis()->GetBinLowEdge(1), proj2->GetXaxis()->GetBinLowEdge(1), division->GetXaxis()->GetBinUpEdge(division->GetNbinsX()), proj2->GetXaxis()->GetBinUpEdge(proj2->GetNbinsX()));
29771dc8	180	//printf("doing axis: %s, y axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsY(), proj2->GetNbinsY(), division->GetYaxis()->GetBinLowEdge(1), proj2->GetYaxis()->GetBinLowEdge(1), division->GetYaxis()->GetBinUpEdge(division->GetNbinsY()), proj2->GetYaxis()->GetBinUpEdge(proj2->GetNbinsY()));
9e952c39	181	division->Divide(proj, proj2, 1, 1, "B");
29771dc8	182	division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle()));
0ab29cfa	183
	184	if (putErrors)
	185	{
	186	division->Sumw2();
	187	if (division->GetDimension() == 1)
	188	{
	189	Int_t nBins = division->GetNbinsX();
29771dc8	190	for (Int_t i = 1; i <= nBins; ++i)
0ab29cfa	191	if (proj2->GetBinContent(i) != 0)
	192	division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i));
	193	}
	194	else if (division->GetDimension() == 2)
	195	{
	196	Int_t nBinsX = division->GetNbinsX();
	197	Int_t nBinsY = division->GetNbinsY();
29771dc8	198	for (Int_t i = 1; i <= nBinsX; ++i)
29771dc8	199	for (Int_t j = 1; j <= nBinsY; ++j)
0ab29cfa	200	if (proj2->GetBinContent(i, j) != 0)
	201	division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
	202	}
	203	}
29771dc8	204
	205	if (save)
	206	{
	207	proj->Write();
	208	proj2->Write();
	209	division->Write();
	210	}
1afae8ff	211	}
1afae8ff	212
92d2d8ad	213	//____________________________________________________________________
25db2d85	214	const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
92d2d8ad	215	{
	216	// returns the title of the axis given in axis (x, y, z)
	217
	218	if (axis == 'x')
	219	return hist->GetXaxis()->GetTitle();
	220	else if (axis == 'y')
	221	return hist->GetYaxis()->GetTitle();
	222	else if (axis == 'z')
	223	return hist->GetZaxis()->GetTitle();
	224
	225	return 0;
	226	}