[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetEmbeddingFromPYTHIATask.cxx

// $Id: AliJetEmbeddingFromPYTHIATask.cxx $
//
// Jet embedding from PYTHIA task.
//
// Author: S.Aiola, C.Loizides

#include "AliJetEmbeddingFromPYTHIATask.h"

#include <TFile.h>
#include <TMath.h>
#include <TString.h>
#include <TRandom.h>
#include <TParameter.h>
#include <TSystem.h>
#include <TH1F.h>

#include "AliVEvent.h"
#include "AliLog.h"

ClassImp(AliJetEmbeddingFromPYTHIATask)

//________________________________________________________________________
AliJetEmbeddingFromPYTHIATask::AliJetEmbeddingFromPYTHIATask() : 
  AliJetEmbeddingFromAODTask("AliJetEmbeddingFromPYTHIATask"),
  fPYTHIAPath(),
  fPtHardBinScaling(),
  fLHC11hAnchorRun(kTRUE),
  fAnchorRun(-1),
  fCurrentPtHardBin(-1),
  fPtHardBinParam(0),
  fHistPtHardBins(0)
{
  // Default constructor.
  SetSuffix("PYTHIAEmbedding");
  fTotalFiles = 2000;
  fRandomAccess = kTRUE;
  SetMC(kTRUE);
}

//________________________________________________________________________
AliJetEmbeddingFromPYTHIATask::AliJetEmbeddingFromPYTHIATask(const char *name, Bool_t drawqa) : 
  AliJetEmbeddingFromAODTask(name, drawqa),
  fPYTHIAPath("/alice/sim/2012/LHC12a15e"),
  fPtHardBinScaling(),
  fLHC11hAnchorRun(kTRUE),
  fAnchorRun(-1),
  fCurrentPtHardBin(-1),
  fPtHardBinParam(0),
  fHistPtHardBins(0)
{
  // Standard constructor.
  SetSuffix("PYTHIAEmbedding");
  fTotalFiles = 2000;
  fRandomAccess = kTRUE;
  SetMC(kTRUE);
}

//________________________________________________________________________
AliJetEmbeddingFromPYTHIATask::~AliJetEmbeddingFromPYTHIATask()
{
  // Destructor
}

//________________________________________________________________________
void AliJetEmbeddingFromPYTHIATask::UserCreateOutputObjects()
{
  if (!fQAhistos)
    return;

  AliJetModelBaseTask::UserCreateOutputObjects();

  fHistPtHardBins = new TH1F("fHistPtHardBins", "fHistPtHardBins", 11, 0, 11);
  fHistPtHardBins->GetXaxis()->SetTitle("p_{T} hard bin");
  fHistPtHardBins->GetYaxis()->SetTitle("total events");
  fOutput->Add(fHistPtHardBins);

  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};

  for (Int_t i = 1; i < 12; i++) 
    fHistPtHardBins->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));

  PostData(1, fOutput);
}

//________________________________________________________________________
Bool_t AliJetEmbeddingFromPYTHIATask::ExecOnce() 
{
  if (fPtHardBinScaling.GetSize() > 0) {
    Double_t sum = 0;
    for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++) 
      sum += fPtHardBinScaling[i];
    
    if (sum == 0) {
      AliWarning("No hard pt bin scaling!");
      sum = fPtHardBinScaling.GetSize();
    }
    
    for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++) 
      fPtHardBinScaling[i] /= sum;
  }

  fPtHardBinParam = static_cast<TParameter<int>*>(InputEvent()->FindListObject("PYTHIAPtHardBin"));
  if (!fPtHardBinParam) {
    fPtHardBinParam = new TParameter<int>("PYTHIAPtHardBin", 0);
    AliDebug(2,"Adding pt hard bin param object to the event list...");
    InputEvent()->AddObject(fPtHardBinParam);
  }

  return AliJetEmbeddingFromAODTask::ExecOnce();
}

//________________________________________________________________________
Bool_t AliJetEmbeddingFromPYTHIATask::GetNextEntry()
{
  Int_t newPtHard = GetRandomPtHardBin();

  new (fPtHardBinParam) TParameter<int>("PYTHIAPtHardBin", newPtHard);

  if (fHistPtHardBins)
    fHistPtHardBins->SetBinContent(newPtHard+1, fHistPtHardBins->GetBinContent(newPtHard+1)+1);

  if (newPtHard != fCurrentPtHardBin) {
    fCurrentPtHardBin = newPtHard;
    if (!OpenNextFile()) return kFALSE;
  }

  return AliJetEmbeddingFromAODTask::GetNextEntry();
}

//________________________________________________________________________
Int_t AliJetEmbeddingFromPYTHIATask::GetRandomPtHardBin() 
{
  static Int_t order[20]={-1};

  if (order[0] == -1)
    TMath::Sort(fPtHardBinScaling.GetSize(), fPtHardBinScaling.GetArray(), order);

  Double_t rnd = gRandom->Rndm();
  Double_t sum = 0;
  Int_t ptHard = -1;
  for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++) {
    sum += fPtHardBinScaling[order[i]];
    if (sum >= rnd) {
      ptHard = order[i];
      break;
    }
  }

  return ptHard;
}

//________________________________________________________________________
Bool_t AliJetEmbeddingFromPYTHIATask::UserNotify()
{
  if (!fLHC11hAnchorRun)
    return kTRUE;
  
  Int_t runNumber = InputEvent()->GetRunNumber();

  Int_t semiGoodRunList[28] = {169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088, 
			       170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195, 
			       170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269, 
			       170270, 170306, 170308, 170309};

  fAnchorRun = 169838; // Assume it is a good run

  for (Int_t i = 0; i < 28; i++) {
    if (runNumber == semiGoodRunList[i]) { // If it is semi good, change the anchor run
      fAnchorRun = 170040;
      break;
    }
  }

  return kTRUE;
}

//________________________________________________________________________
TString AliJetEmbeddingFromPYTHIATask::GetNextFileName() 
{
  fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*(fTotalFiles-1))+1;

  TString fileName(fPYTHIAPath);
  if (!fileName.EndsWith("/"))
    fileName += "/";
  
  fileName += fAnchorRun;
  fileName += "/";
  fileName += fCurrentPtHardBin;
  fileName += "/";
  if (fCurrentAODFileID < 10)
    fileName += "00";
  else if (fCurrentAODFileID < 100)
    fileName += "0";
  fileName += fCurrentAODFileID;

  if (!gSystem->AccessPathName(Form("%s/AliAOD.root", fileName.Data())))
    fileName += "/AliAOD.root";
  else if (!gSystem->AccessPathName(Form("%s/aod_archive.zip", fileName.Data())))
    fileName += "/aod_archive.zip#AliAOD.root";
  else
    fileName += "/root_archive.zip#AliAOD.root";

  return fileName;
}
Commit	Line	Data
a6e51e61	1	// $Id: AliJetEmbeddingFromPYTHIATask.cxx $
	2	//
	3	// Jet embedding from PYTHIA task.
	4	//
	5	// Author: S.Aiola, C.Loizides
	6
	7	#include "AliJetEmbeddingFromPYTHIATask.h"
	8
	9	#include <TFile.h>
	10	#include <TMath.h>
	11	#include <TString.h>
	12	#include <TRandom.h>
ca5c29fa	13	#include <TParameter.h>
	14	#include <TSystem.h>
	15	#include <TH1F.h>
a6e51e61	16
	17	#include "AliVEvent.h"
	18	#include "AliLog.h"
	19
	20	ClassImp(AliJetEmbeddingFromPYTHIATask)
	21
	22	//________________________________________________________________________
	23	AliJetEmbeddingFromPYTHIATask::AliJetEmbeddingFromPYTHIATask() :
ca5c29fa	24	AliJetEmbeddingFromAODTask("AliJetEmbeddingFromPYTHIATask"),
a6e51e61	25	fPYTHIAPath(),
	26	fPtHardBinScaling(),
	27	fLHC11hAnchorRun(kTRUE),
	28	fAnchorRun(-1),
ca5c29fa	29	fCurrentPtHardBin(-1),
	30	fPtHardBinParam(0),
	31	fHistPtHardBins(0)
a6e51e61	32	{
	33	// Default constructor.
	34	SetSuffix("PYTHIAEmbedding");
	35	fTotalFiles = 2000;
7030f36f	36	fRandomAccess = kTRUE;
5ce8ae64	37	SetMC(kTRUE);
a6e51e61	38	}
	39
	40	//________________________________________________________________________
ca5c29fa	41	AliJetEmbeddingFromPYTHIATask::AliJetEmbeddingFromPYTHIATask(const char *name, Bool_t drawqa) :
ca5c29fa	42	AliJetEmbeddingFromAODTask(name, drawqa),
a6e51e61	43	fPYTHIAPath("/alice/sim/2012/LHC12a15e"),
	44	fPtHardBinScaling(),
	45	fLHC11hAnchorRun(kTRUE),
	46	fAnchorRun(-1),
ca5c29fa	47	fCurrentPtHardBin(-1),
	48	fPtHardBinParam(0),
	49	fHistPtHardBins(0)
a6e51e61	50	{
	51	// Standard constructor.
	52	SetSuffix("PYTHIAEmbedding");
	53	fTotalFiles = 2000;
7030f36f	54	fRandomAccess = kTRUE;
5ce8ae64	55	SetMC(kTRUE);
a6e51e61	56	}
	57
	58	//________________________________________________________________________
	59	AliJetEmbeddingFromPYTHIATask::~AliJetEmbeddingFromPYTHIATask()
	60	{
	61	// Destructor
	62	}
	63
ca5c29fa	64	//________________________________________________________________________
	65	void AliJetEmbeddingFromPYTHIATask::UserCreateOutputObjects()
	66	{
	67	if (!fQAhistos)
	68	return;
	69
	70	AliJetModelBaseTask::UserCreateOutputObjects();
	71
	72	fHistPtHardBins = new TH1F("fHistPtHardBins", "fHistPtHardBins", 11, 0, 11);
	73	fHistPtHardBins->GetXaxis()->SetTitle("p_{T} hard bin");
	74	fHistPtHardBins->GetYaxis()->SetTitle("total events");
	75	fOutput->Add(fHistPtHardBins);
	76
	77	const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
	78	const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
	79
	80	for (Int_t i = 1; i < 12; i++)
	81	fHistPtHardBins->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
	82
	83	PostData(1, fOutput);
	84	}
	85
a6e51e61	86	//________________________________________________________________________
	87	Bool_t AliJetEmbeddingFromPYTHIATask::ExecOnce()
	88	{
	89	if (fPtHardBinScaling.GetSize() > 0) {
	90	Double_t sum = 0;
	91	for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++)
	92	sum += fPtHardBinScaling[i];
	93
	94	if (sum == 0) {
	95	AliWarning("No hard pt bin scaling!");
	96	sum = fPtHardBinScaling.GetSize();
	97	}
	98
	99	for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++)
	100	fPtHardBinScaling[i] /= sum;
	101	}
	102
ca5c29fa	103	fPtHardBinParam = static_cast<TParameter<int>*>(InputEvent()->FindListObject("PYTHIAPtHardBin"));
	104	if (!fPtHardBinParam) {
	105	fPtHardBinParam = new TParameter<int>("PYTHIAPtHardBin", 0);
	106	AliDebug(2,"Adding pt hard bin param object to the event list...");
	107	InputEvent()->AddObject(fPtHardBinParam);
	108	}
	109
a6e51e61	110	return AliJetEmbeddingFromAODTask::ExecOnce();
	111	}
	112
	113	//________________________________________________________________________
63fac07f	114	Bool_t AliJetEmbeddingFromPYTHIATask::GetNextEntry()
a6e51e61	115	{
	116	Int_t newPtHard = GetRandomPtHardBin();
	117
63fac07f	118	new (fPtHardBinParam) TParameter<int>("PYTHIAPtHardBin", newPtHard);
ca5c29fa	119
ca5c29fa	120	if (fHistPtHardBins)
63fac07f	121	fHistPtHardBins->SetBinContent(newPtHard+1, fHistPtHardBins->GetBinContent(newPtHard+1)+1);
ca5c29fa	122
a6e51e61	123	if (newPtHard != fCurrentPtHardBin) {
	124	fCurrentPtHardBin = newPtHard;
	125	if (!OpenNextFile()) return kFALSE;
	126	}
	127
	128	return AliJetEmbeddingFromAODTask::GetNextEntry();
	129	}
	130
	131	//________________________________________________________________________
	132	Int_t AliJetEmbeddingFromPYTHIATask::GetRandomPtHardBin()
	133	{
	134	static Int_t order[20]={-1};
	135
	136	if (order[0] == -1)
	137	TMath::Sort(fPtHardBinScaling.GetSize(), fPtHardBinScaling.GetArray(), order);
	138
	139	Double_t rnd = gRandom->Rndm();
	140	Double_t sum = 0;
	141	Int_t ptHard = -1;
	142	for (Int_t i = 0; i < fPtHardBinScaling.GetSize(); i++) {
	143	sum += fPtHardBinScaling[order[i]];
	144	if (sum >= rnd) {
	145	ptHard = order[i];
	146	break;
	147	}
	148	}
	149
	150	return ptHard;
	151	}
	152
	153	//________________________________________________________________________
	154	Bool_t AliJetEmbeddingFromPYTHIATask::UserNotify()
	155	{
	156	if (!fLHC11hAnchorRun)
	157	return kTRUE;
	158
	159	Int_t runNumber = InputEvent()->GetRunNumber();
	160
	161	Int_t semiGoodRunList[28] = {169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088,
	162	170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195,
	163	170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269,
	164	170270, 170306, 170308, 170309};
	165
	166	fAnchorRun = 169838; // Assume it is a good run
	167
	168	for (Int_t i = 0; i < 28; i++) {
	169	if (runNumber == semiGoodRunList[i]) { // If it is semi good, change the anchor run
	170	fAnchorRun = 170040;
	171	break;
	172	}
	173	}
	174
	175	return kTRUE;
	176	}
	177
	178	//________________________________________________________________________
	179	TString AliJetEmbeddingFromPYTHIATask::GetNextFileName()
	180	{
	181	fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*(fTotalFiles-1))+1;
	182
	183	TString fileName(fPYTHIAPath);
	184	if (!fileName.EndsWith("/"))
	185	fileName += "/";
	186
187	fileName += fAnchorRun;
188	fileName += "/";
189	fileName += fCurrentPtHardBin;
190	fileName += "/";
191	if (fCurrentAODFileID < 10)
192	fileName += "00";
193	else if (fCurrentAODFileID < 100)
194	fileName += "0";
195	fileName += fCurrentAODFileID;
ca5c29fa	196
63fac07f	197	if (!gSystem->AccessPathName(Form("%s/AliAOD.root", fileName.Data())))
ca5c29fa	198	fileName += "/AliAOD.root";
63fac07f	199	else if (!gSystem->AccessPathName(Form("%s/aod_archive.zip", fileName.Data())))
ca5c29fa	200	fileName += "/aod_archive.zip#AliAOD.root";
	201	else
	202	fileName += "/root_archive.zip#AliAOD.root";
a6e51e61	203
	204	return fileName;
	205	}