[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, Trigger trigger)
{
  // 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();

  switch (trigger)
  {
    case kMB1:
    {
      if (triggerMask&32 || ((triggerMask&1) || (triggerMask&2)))
        return kTRUE;
      break;
    }
    case kMB2:
    {
      if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2)))
        return kTRUE;
      break;
    }
  }

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