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

  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(proj2);
  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
	59	return kTRUE;
	60	}
	61
	62	//____________________________________________________________________
7584d357	63	Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
04a7657f	64	{
04a7657f	65	//
25db2d85	66	// this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
	67	// shall be counted as a primary particle
	68	//
04a7657f	69	// This function or a equivalent should be available in some common place of AliRoot
	70	//
	71
	72	// if the particle has a daughter primary, we do not want to count it
	73	if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
	74	{
7584d357	75	if (adebug)
25db2d85	76	printf("Dropping particle because it has a daughter among the primaries.\n");
04a7657f	77	return kFALSE;
	78	}
	79
	80	Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
	81
	82	// skip quarks and gluon
	83	if (pdgCode <= 10 \|\| pdgCode == 21)
	84	{
7584d357	85	if (adebug)
25db2d85	86	printf("Dropping particle because it is a quark or gluon.\n");
04a7657f	87	return kFALSE;
	88	}
	89
	90	if (strcmp(aParticle->GetName(),"XXX") == 0)
	91	{
7584d357	92	if (adebug)
25db2d85	93	printf("WARNING: There is a particle named XXX.\n");
04a7657f	94	return kFALSE;
	95	}
	96
	97	TParticlePDG* pdgPart = aParticle->GetPDG();
	98
	99	if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
	100	{
7584d357	101	if (adebug)
25db2d85	102	printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode);
04a7657f	103	return kFALSE;
	104	}
	105
	106	if (pdgPart->Charge() == 0)
	107	{
7584d357	108	if (adebug)
25db2d85	109	printf("Dropping particle because it is not charged.\n");
04a7657f	110	return kFALSE;
	111	}
	112
	113	return kTRUE;
	114	}
847489f7	115
847489f7	116	//____________________________________________________________________
29771dc8	117	void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save)
847489f7	118	{
	119	// create projections of 3d hists to all 2d combinations
	120	// the histograms are not returned, just use them from memory or use this to create them in a file
	121
	122	TH1* proj = hist->Project3D("yx");
	123	proj->SetXTitle(hist->GetXaxis()->GetTitle());
	124	proj->SetYTitle(hist->GetYaxis()->GetTitle());
29771dc8	125	if (save)
29771dc8	126	proj->Write();
847489f7	127
	128	proj = hist->Project3D("zx");
	129	proj->SetXTitle(hist->GetXaxis()->GetTitle());
	130	proj->SetYTitle(hist->GetZaxis()->GetTitle());
29771dc8	131	if (save)
29771dc8	132	proj->Write();
847489f7	133
	134	proj = hist->Project3D("zy");
	135	proj->SetXTitle(hist->GetYaxis()->GetTitle());
	136	proj->SetYTitle(hist->GetZaxis()->GetTitle());
29771dc8	137	if (save)
29771dc8	138	proj->Write();
847489f7	139	}
1afae8ff	140
1afae8ff	141	//____________________________________________________________________
29771dc8	142	void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save)
1afae8ff	143	{
	144	// create projections of the 3d hists divides them
	145	// axis decides to which plane, if axis is 0 to all planes
	146	// the histograms are not returned, just use them from memory or use this to create them in a file
	147
	148	if (axis == 0)
	149	{
29771dc8	150	CreateDividedProjections(hist, hist2, "yx", putErrors, save);
	151	CreateDividedProjections(hist, hist2, "zx", putErrors, save);
	152	CreateDividedProjections(hist, hist2, "zy", putErrors, save);
1afae8ff	153
	154	return;
	155	}
	156
	157	TH1* proj = hist->Project3D(axis);
0ab29cfa	158
	159	if (strlen(axis) == 2)
	160	{
	161	proj->SetYTitle(GetAxisTitle(hist, axis[0]));
	162	proj->SetXTitle(GetAxisTitle(hist, axis[1]));
	163	}
	164	else if (strlen(axis) == 1)
	165	proj->SetXTitle(GetAxisTitle(hist, axis[0]));
1afae8ff	166
1afae8ff	167	TH1* proj2 = hist2->Project3D(axis);
0ab29cfa	168	if (strlen(axis) == 2)
	169	{
	170	proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
	171	proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
	172	}
	173	else if (strlen(axis) == 1)
	174	proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));
1afae8ff	175
1afae8ff	176	TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
29771dc8	177	//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	178	//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()));
1afae8ff	179	division->Divide(proj2);
29771dc8	180	division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle()));
0ab29cfa	181
	182	if (putErrors)
	183	{
	184	division->Sumw2();
	185	if (division->GetDimension() == 1)
	186	{
	187	Int_t nBins = division->GetNbinsX();
29771dc8	188	for (Int_t i = 1; i <= nBins; ++i)
0ab29cfa	189	if (proj2->GetBinContent(i) != 0)
	190	division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i));
	191	}
	192	else if (division->GetDimension() == 2)
	193	{
	194	Int_t nBinsX = division->GetNbinsX();
	195	Int_t nBinsY = division->GetNbinsY();
29771dc8	196	for (Int_t i = 1; i <= nBinsX; ++i)
29771dc8	197	for (Int_t j = 1; j <= nBinsY; ++j)
0ab29cfa	198	if (proj2->GetBinContent(i, j) != 0)
	199	division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
	200	}
	201	}
29771dc8	202
	203	if (save)
	204	{
	205	proj->Write();
	206	proj2->Write();
	207	division->Write();
	208	}
1afae8ff	209	}
1afae8ff	210
92d2d8ad	211	//____________________________________________________________________
25db2d85	212	const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
92d2d8ad	213	{
	214	// returns the title of the axis given in axis (x, y, z)
	215
	216	if (axis == 'x')
	217	return hist->GetXaxis()->GetTitle();
	218	else if (axis == 'y')
	219	return hist->GetYaxis()->GetTitle();
	220	else if (axis == 'z')
	221	return hist->GetZaxis()->GetTitle();
	222
	223	return 0;
	224	}