[u/mrichter/AliRoot.git] / HLT / BASE / util / AliHLTMCGeneratorComponent.cxx

// -*- Mode: C++ -*-
// $Id: AliHLTMCGeneratorComponent.cxx  $
/**************************************************************************
 * This file is property of and copyright by the ALICE HLT Project        * 
 * ALICE Experiment at CERN, All rights reserved.                         *
 *                                                                        *
 * Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de>        *
 *                  for The ALICE HLT Project.                            *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/** @file   AliHLTMCGeneratorComponent.cxx
    @author Jochen Thaeder
    @date   
    @brief  Component for generating MC events
    @note   The class is used in Offline (AliRoot) context
*/

#include "AliHLTMCGeneratorComponent.h"

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TH1F.h>
#include <TStopwatch.h>
#include <TDatime.h>
#include <TRandom.h>
#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TArrayI.h>

#include "AliGenerator.h"
#include "AliPDG.h"
#include "AliRunLoader.h"
#include "AliRun.h"
#include "AliStack.h"
#include "AliHeader.h"
#include "AliGenPythia.h"
#include "AliPythia.h"
#endif

#include "TPythia6Calls.h"
#include "TPythia6.h"

/** ROOT macro for the implementation of ROOT specific class methods */
ClassImp(AliHLTMCGeneratorComponent)

/*
 * ---------------------------------------------------------------------------------
 *                            Constructor / Destructor
 * ---------------------------------------------------------------------------------
 */

// #################################################################################
AliHLTMCGeneratorComponent::AliHLTMCGeneratorComponent() :  
  AliHLTDataSource(),
  fpMC(NULL),
  fpHLTMC(NULL),
  fNEvents(0),
  fCurrentEvent(0),
  fEventNumber(0),
  fRunNumber(0),
  fRunLoader(NULL),
  fGenerator(NULL),
  fComment(""),
  fSeed(0),
  fRunType(kPythia6Jets104_125),
  fEcms(14000.),
  fJetEtaMax(0.2), 
  fJetEtaMin(-0.2),
  fJetEtMax(1000.),
  fJetEtMin(10.0),
  fJetConeRadius(0.4),
  fPtHardMin(10.0),
  fPtHardMax(-1.0),
  fQuenching(0),                   
  fQhat(20.),
  fApplyParticleCuts(kFALSE) {
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
}

// #################################################################################
AliHLTMCGeneratorComponent::~AliHLTMCGeneratorComponent() {
  // see header file for class documentation

  // file list and file name list are owner of their objects and
  // delete all the objects
}

/*
 * ---------------------------------------------------------------------------------
 *                              Initialize static const
 * ---------------------------------------------------------------------------------
 */

//##################################################################################
const Char_t *AliHLTMCGeneratorComponent::fgkPprRunName[] = {
  "kPythia6Jets20_24",   "kPythia6Jets24_29",   "kPythia6Jets29_35",
  "kPythia6Jets35_42",   "kPythia6Jets42_50",   "kPythia6Jets50_60",
  "kPythia6Jets60_72",   "kPythia6Jets72_86",   "kPythia6Jets86_104",
  "kPythia6Jets104_125", "kPythia6Jets125_150", "kPythia6Jets150_180",
  "kPyJetJet", "kPyGammaJetPHOS"
};

/*
 * ---------------------------------------------------------------------------------
 * Public functions to implement AliHLTComponent's interface.
 * These functions are required for the registration process
 * ---------------------------------------------------------------------------------
 */

// #################################################################################
const Char_t* AliHLTMCGeneratorComponent::GetComponentID() {
  // see header file for class documentation
  return "MCGenerator";
}

// #################################################################################
AliHLTComponentDataType AliHLTMCGeneratorComponent::GetOutputDataType() {
  // see header file for class documentation

  return kAliHLTDataTypeMCObject|kAliHLTDataOriginHLT;
}

// #################################################################################
void AliHLTMCGeneratorComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) {
  // see header file for class documentation

  constBase = 20000;
  inputMultiplier = 1.0;
}

// #################################################################################
AliHLTComponent* AliHLTMCGeneratorComponent::Spawn() {
  // see header file for class documentation
  return new AliHLTMCGeneratorComponent;
}

/*
 * ---------------------------------------------------------------------------------
 * Protected functions to implement AliHLTComponent's interface.
 * These functions provide initialization as well as the actual processing
 * capabilities of the component. 
 * ---------------------------------------------------------------------------------
 */

// #################################################################################
Int_t AliHLTMCGeneratorComponent::DoInit(int argc, const char** argv) {
  // see header file for class documentation

  Int_t iResult = 0;
  Int_t bMissingParam=0;
  
  TString argument="";

  // -- Loop over all arguments
  for ( Int_t iter = 0; iter<argc && iResult>=0; iter++) {
    argument=argv[iter];
    if (argument.IsNull()) continue;
    
    // -- seed
    if ( !argument.CompareTo("-seed") ) {
      if ((bMissingParam=(++iter>=argc))) break;
      
      TString parameter(argv[iter]);
      parameter.Remove(TString::kLeading, ' ');

      if ( parameter.IsDigit() ) {
	fSeed = parameter.Atoi();
      }
      
      else {
	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
	iResult=-EINVAL;
      }
    } 

    // -- runtype
    else if ( !argument.CompareTo("-runtype") ) {
      if ((bMissingParam=(++iter>=argc))) break;
      
      TString parameter(argv[iter]);
      parameter.Remove(TString::kLeading, ' ');

      Bool_t found = kFALSE;

      if ( parameter.IsDigit() ) {
	Int_t idx = parameter.Atoi();
	if ( idx < kRunMax ) {
	  fRunType = (PprRun_t) idx;
	  found = kTRUE;
	}
      }
      else {
	for ( Int_t iRun = 0; iRun < kRunMax && !found; iRun++ ) {
	  if ( ! parameter.CompareTo(fgkPprRunName[iRun]) ) {
	    fRunType = (PprRun_t) iRun;  
	    found = kTRUE;
	  }
	}      
      }
      if ( !found ) {
	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
	iResult=-EINVAL;
      }
    }

    // -- nevents
    else if ( !argument.CompareTo("-nevents") ) {
      if ((bMissingParam=(++iter>=argc))) break;
      
      TString parameter(argv[iter]);
      parameter.Remove(TString::kLeading, ' ');

      if ( parameter.IsDigit() ) {
	fNEvents = parameter.Atoi();
      }
      
      else {
	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
	iResult=-EINVAL;
      }
    } 

    // -- coneRadius
    else if ( !argument.CompareTo("-coneRadius") ) {
      if ((bMissingParam=(++iter>=argc))) break;
      
      TString parameter(argv[iter]);
      parameter.Remove(TString::kLeading, ' ');
      
      if ( parameter.IsFloat() ) {
	fJetConeRadius = parameter.Atof();
      }
      else {
	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
	iResult=-EINVAL;
      }
    } 

    // -- jetCutMinEt
    else if ( !argument.CompareTo("-jetCutMinEt") ) {
      if ((bMissingParam=(++iter>=argc))) break;
      
      TString parameter(argv[iter]);
      parameter.Remove(TString::kLeading, ' ');
      
      if ( parameter.IsFloat() ) {
	fJetEtMin = parameter.Atof();
      }
      else {
	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
	iResult=-EINVAL;
      }
    } 

    // -- applyParticleCuts
    else if ( !argument.CompareTo("-applyParticleCuts") ) {
      fApplyParticleCuts = kTRUE;
    } 

    // -- Argument not known
    else {
      HLTError("Unknown argument %s.", argument.Data());
      iResult = -EINVAL;
    }
    
  } // for ( Int iter = 0; iter<argc && iResult>=0; iter++) {
  
  // -- Check if parameter is missing
  if ( bMissingParam ) {
    HLTError("Missing parameter for argument %s.", argument.Data());
     iResult=-EINVAL;
  }

  //  -- Check if seed and nevents has been specified
  if ( !fSeed || !fNEvents ) {
    HLTError("-seed and -nevents need to be specified.");
    iResult=-EINVAL;
  }

  // -- Set Random Number seed
  gRandom->SetSeed( fSeed );

  // -- Run loader 
  // ---------------
  AliPDG::AddParticlesToPdgDataBase();

  fRunLoader = AliRunLoader::Open("galice.root","FASTRUN","recreate");
  if (fRunLoader == NULL ) {
    //    gAlice->Fatal("AliHLTMCGeneratorComponent","Can not instatiate the Run Loader");
    HLTFatal("Can not instatiate the Run Loader");
    return -EINPROGRESS;
  }

  fRunLoader->SetCompressionLevel(2);
  fRunLoader->SetNumberOfEventsPerFile(fNEvents);
  fRunLoader->LoadKinematics("RECREATE");
  fRunLoader->MakeTree("E");
  gAlice->SetRunLoader(fRunLoader); 

  // -- Create stack
  fRunLoader->MakeStack();
 
  // -- Create and Initialize Generator
  fGenerator = GeneratorFactory();
  fGenerator->SetStack(fRunLoader->Stack());
  fGenerator->Init();

  HLTInfo("MC Generator setup with : %s", fComment.Data() );

  return iResult;
}

// #################################################################################
Int_t AliHLTMCGeneratorComponent::DoDeinit() {
  // see header file for class documentation

  if ( fpMC ) 
    delete fpMC;
  fpMC = NULL;

  if ( fpHLTMC ) 
    delete fpHLTMC;
  fpHLTMC = NULL;

  return 0;
}

// #################################################################################
Int_t AliHLTMCGeneratorComponent::GetEvent( const AliHLTComponentEventData& /*evtData*/,
						AliHLTComponentTriggerData& /*trigData*/,
						AliHLTUInt8_t* /*outputPtr*/, 
						AliHLTUInt32_t& size,
						vector<AliHLTComponentBlockData>& /*outputBlocks*/ ) {
  // see header file for class documentation

  if ( !IsDataEvent() ) 
    return 0;

  Int_t iResult=0;
  size=0;

  AliStack*  stack  = fRunLoader->Stack();
  AliHeader* header = fRunLoader->GetHeader();
  
  // -- Initialize event
  header->Reset(0,fCurrentEvent);

  // -- Reset stack 
  stack->Reset();

  // -- Generate event
  fGenerator->Generate();
    
  // -- Finish event
  header->SetNprimary(stack->GetNprimary());
  header->SetNtrack(stack->GetNtrack());  

  // -- Finish Event	
  stack->FinishEvent();
  header->SetStack(stack);

  // -- Create HLT MC Event
  if ( fpHLTMC ) 
    delete fpHLTMC;
  fpHLTMC = new AliHLTMCEvent( fApplyParticleCuts);

  // -- Fill HLT MC event
  iResult = fpHLTMC->FillMCEvent( stack, header );

  // -- Send HLT MC event
  if ( fpHLTMC && !iResult )
    PushBack( fpHLTMC, kAliHLTDataTypeMCObject|kAliHLTDataOriginHLT, 0 ); 

  fCurrentEvent++;

  return iResult;
}

// #################################################################################
AliGenerator* AliHLTMCGeneratorComponent::GeneratorFactory() {
  
  // -- Configuration for hard QCD processes generation with PYTHIA 
  // ----------------------------------------------------------------
  AliGenPythia * gener = new AliGenPythia(-1);

  gener->SetEnergyCMS(fEcms); // Centre of mass energy
  gener->SetGluonRadiation(1,1);
  gener->SetForceDecay(kAll); //  Decay type

  // -- Final state kinematic cuts
  gener->SetJetEtaRange(fJetEtaMin, fJetEtaMax);
  gener->SetJetPhiRange(0., 360.);
  gener->SetJetEtRange(fJetEtMin, fJetEtMax);
  // gener->SetMomentumRange(0,99999999);

  // -- Structure function
  gener->SetStrucFunc(kCTEQ4L);

  /************************************************************************
   * void AliGenPythia::SetPycellParameters()
   *  Float_t etamax   -- fPycellEtaMax     -- PARU(51) -- +/- eta max of detector
   *  Int_t neta       -- fPycellNEta       -- MSTU(51) -- N cells in eta
   *  Int_t nphi       -- fPycellNPhi       -- MSTU(52) -- N cells in phi
   *  Float_t thresh   -- fPycellThreshold  -- PARU(58) -- cells with Et below treshold are discared
   *  Float_t etseed   -- fPycellEtSeed     -- PARU(52) -- Et threshold for seed cells
   *  Float_t minet    -- fPycellMinEtJet   -- PARU(53) -- Min Et for jets
   *  Float_t r        -- fPycellMaxRadius  -- PARU(54) -- cone radius
   ************************************************************************/ 
  gener->SetPycellParameters(fJetEtaMax, 274, 432, 0., 4., fJetEtMin, fJetConeRadius);

  // gener->SetPtKick(5.);      // set the intrinsic kt to 5 GeV/c
  
  // -- Jet Quenching
  gener->SetQuench(fQuenching);
  if( fQuenching == 1){
    Float_t k = 6e5*(fQhat/1.7) ; // qhat=1.7, k = 6e5, default  value 
    AliPythia::Instance()->InitQuenching(0.,0.1,k,0,0.95,6);
  }

  switch (fRunType) {
  case kPythia6Jets20_24:
    {
      fComment = fComment.Append(":Pythia jets 20-24 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(20., 24.);// Pt transfer of the hard scattering
      fPtHardMin=20.;
    }
    break;
  case kPythia6Jets24_29:
    {
      fComment = fComment.Append(":Pythia jets 24-29 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(24., 29.);// Pt transfer of the hard scattering
      fPtHardMin=24.;
    }
    break;
  case kPythia6Jets29_35:
    {
      fComment = fComment.Append(":Pythia jets 29-35 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(29., 35.);// Pt transfer of the hard scattering
      fPtHardMin=29.;
    }
    break;
  case kPythia6Jets35_42:
    {
      fComment = fComment.Append(":Pythia jets 35-42 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(35., 42.);// Pt transfer of the hard scattering
      fPtHardMin=35.;
    }
    break;
  case kPythia6Jets42_50:
    {
      fComment = fComment.Append(":Pythia jets 42-50 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(42., 50.);// Pt transfer of the hard scattering
      fPtHardMin=42.;
    }
    break;
  case kPythia6Jets50_60:
    {
      fComment = fComment.Append(":Pythia jets 50-60 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(50., 60.);// Pt transfer of the hard scattering
      fPtHardMin=50.;
    }
    break;
  case kPythia6Jets60_72:
    {
      fComment = fComment.Append(":Pythia jets 60-72 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(60., 72.);// Pt transfer of the hard scattering
      fPtHardMin=60.;
    }
    break;
  case kPythia6Jets72_86:
    {
      fComment = fComment.Append(":Pythia jets 72-86 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(72., 86.);// Pt transfer of the hard scattering
      fPtHardMin=72.;
    }
    break;
  case kPythia6Jets86_104:
    {
      fComment = fComment.Append(":Pythia jets 86-104 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(86., 104.);// Pt transfer of the hard scattering
      fPtHardMin=86.;
    }
    break;
  case kPythia6Jets104_125:
    {
      fComment = fComment.Append(":Pythia jets 104-125 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(104., 125.);// Pt transfer of the hard scattering
      fPtHardMin=104.;
    }
    break;
  case kPythia6Jets125_150:
    {
      fComment = fComment.Append(":Pythia jets 125-150 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(125., 150.);// Pt transfer of the hard scattering
      fPtHardMin=125.;
    }
    break;
  case kPythia6Jets150_180:
    {
      fComment = fComment.Append(":Pythia jets 150-180 GeV @ 14 TeV");
      gener->SetProcess(kPyJets);//        Process type
      gener->SetPtHard(150., 180.);// Pt transfer of the hard scattering
      fPtHardMin=150.;
    }
    break;
  case kPyJetJet:
    {
      fComment = fComment.Append(" Jet-jet @ 14 TeV");
      gener->SetProcess(kPyJets); 
      gener->SetPtHard(fPtHardMin,fPtHardMax);
      gener->SetEventListRange(0,1);  // XXXX
    }
    break;
  case kPyGammaJetPHOS:
    {
      fComment = fComment.Append(" pp->jet + Gamma over PHOS @ 14 TeV");
      gener->SetProcess(kPyDirectGamma);
      gener->SetPtHard(fPtHardMin,fPtHardMax);
      //gener->SetYHard(-1.0,1.0);
      gener->SetGammaEtaRange(-0.13,0.13);
      gener->SetGammaPhiRange(218.,322.); //Over 5 modules +-2 deg  
      gener->SetEventListRange(0,1); // XXXX
    }
    break;
    
  default: break;
  }
  
  return gener;
}
Commit	Line	Data
7aac8168	1	// -- Mode: C++ --
	2	// $Id: AliHLTMCGeneratorComponent.cxx $
	3	/**************************************************************************
	4	* This file is property of and copyright by the ALICE HLT Project *
	5	* ALICE Experiment at CERN, All rights reserved. *
	6	* *
	7	* Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de> *
	8	* for The ALICE HLT Project. *
	9	* *
	10	* Permission to use, copy, modify and distribute this software and its *
	11	* documentation strictly for non-commercial purposes is hereby granted *
	12	* without fee, provided that the above copyright notice appears in all *
	13	* copies and that both the copyright notice and this permission notice *
	14	* appear in the supporting documentation. The authors make no claims *
	15	* about the suitability of this software for any purpose. It is *
	16	* provided "as is" without express or implied warranty. *
	17	**************************************************************************/
	18
	19	/** @file AliHLTMCGeneratorComponent.cxx
	20	@author Jochen Thaeder
	21	@date
	22	@brief Component for generating MC events
	23	@note The class is used in Offline (AliRoot) context
	24	*/
	25
	26	#include "AliHLTMCGeneratorComponent.h"
	27
	28	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	29	#include <Riostream.h>
	30	#include <TH1F.h>
	31	#include <TStopwatch.h>
	32	#include <TDatime.h>
	33	#include <TRandom.h>
	34	#include <TDatabasePDG.h>
	35	#include <TParticle.h>
	36	#include <TArrayI.h>
	37
	38	#include "AliGenerator.h"
	39	#include "AliPDG.h"
	40	#include "AliRunLoader.h"
	41	#include "AliRun.h"
	42	#include "AliStack.h"
	43	#include "AliHeader.h"
	44	#include "AliGenPythia.h"
	45	#include "AliPythia.h"
	46	#endif
	47
	48	#include "TPythia6Calls.h"
	49	#include "TPythia6.h"
	50
	51	/** ROOT macro for the implementation of ROOT specific class methods */
	52	ClassImp(AliHLTMCGeneratorComponent)
	53
	54	/*
	55	* ---------------------------------------------------------------------------------
	56	* Constructor / Destructor
	57	* ---------------------------------------------------------------------------------
	58	*/
	59
	60	// #################################################################################
	61	AliHLTMCGeneratorComponent::AliHLTMCGeneratorComponent() :
	62	AliHLTDataSource(),
	63	fpMC(NULL),
	64	fpHLTMC(NULL),
65	fNEvents(0),
66	fCurrentEvent(0),
67	fEventNumber(0),
68	fRunNumber(0),
69	fRunLoader(NULL),
70	fGenerator(NULL),
71	fComment(""),
72	fSeed(0),
73	fRunType(kPythia6Jets104_125),
74	fEcms(14000.),
75	fJetEtaMax(0.2),
76	fJetEtaMin(-0.2),
77	fJetEtMax(1000.),
78	fJetEtMin(10.0),
79	fJetConeRadius(0.4),
80	fPtHardMin(10.0),
81	fPtHardMax(-1.0),
82	fQuenching(0),
396e75ed	83	fQhat(20.),
396e75ed	84	fApplyParticleCuts(kFALSE) {
7aac8168	85	// see header file for class documentation
	86	// or
	87	// refer to README to build package
	88	// or
	89	// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
	90	}
	91
	92	// #################################################################################
	93	AliHLTMCGeneratorComponent::~AliHLTMCGeneratorComponent() {
	94	// see header file for class documentation
	95
	96	// file list and file name list are owner of their objects and
	97	// delete all the objects
	98	}
	99
	100	/*
	101	* ---------------------------------------------------------------------------------
	102	* Initialize static const
	103	* ---------------------------------------------------------------------------------
	104	*/
	105
	106	//##################################################################################
	107	const Char_t *AliHLTMCGeneratorComponent::fgkPprRunName[] = {
	108	"kPythia6Jets20_24", "kPythia6Jets24_29", "kPythia6Jets29_35",
	109	"kPythia6Jets35_42", "kPythia6Jets42_50", "kPythia6Jets50_60",
	110	"kPythia6Jets60_72", "kPythia6Jets72_86", "kPythia6Jets86_104",
	111	"kPythia6Jets104_125", "kPythia6Jets125_150", "kPythia6Jets150_180",
	112	"kPyJetJet", "kPyGammaJetPHOS"
	113	};
	114
	115	/*
	116	* ---------------------------------------------------------------------------------
	117	* Public functions to implement AliHLTComponent's interface.
	118	* These functions are required for the registration process
	119	* ---------------------------------------------------------------------------------
	120	*/
	121
	122	// #################################################################################
	123	const Char_t* AliHLTMCGeneratorComponent::GetComponentID() {
	124	// see header file for class documentation
	125	return "MCGenerator";
	126	}
	127
	128	// #################################################################################
	129	AliHLTComponentDataType AliHLTMCGeneratorComponent::GetOutputDataType() {
	130	// see header file for class documentation
	131
	132	return kAliHLTDataTypeMCObject\|kAliHLTDataOriginHLT;
	133	}
	134
	135	// #################################################################################
86a6fdbc	136	void AliHLTMCGeneratorComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) {
7aac8168	137	// see header file for class documentation
	138
	139	constBase = 20000;
	140	inputMultiplier = 1.0;
	141	}
	142
	143	// #################################################################################
	144	AliHLTComponent* AliHLTMCGeneratorComponent::Spawn() {
	145	// see header file for class documentation
	146	return new AliHLTMCGeneratorComponent;
	147	}
	148
	149	/*
	150	* ---------------------------------------------------------------------------------
	151	* Protected functions to implement AliHLTComponent's interface.
	152	* These functions provide initialization as well as the actual processing
	153	* capabilities of the component.
	154	* ---------------------------------------------------------------------------------
	155	*/
	156
	157	// #################################################################################
	158	Int_t AliHLTMCGeneratorComponent::DoInit(int argc, const char** argv) {
	159	// see header file for class documentation
	160
	161	Int_t iResult = 0;
	162	Int_t bMissingParam=0;
	163
	164	TString argument="";
	165
	166	// -- Loop over all arguments
	167	for ( Int_t iter = 0; iter<argc && iResult>=0; iter++) {
	168	argument=argv[iter];
	169	if (argument.IsNull()) continue;
	170
	171	// -- seed
	172	if ( !argument.CompareTo("-seed") ) {
	173	if ((bMissingParam=(++iter>=argc))) break;
	174
	175	TString parameter(argv[iter]);
	176	parameter.Remove(TString::kLeading, ' ');
	177
	178	if ( parameter.IsDigit() ) {
	179	fSeed = parameter.Atoi();
	180	}
	181
	182	else {
	183	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
	184	iResult=-EINVAL;
	185	}
	186	}
	187
	188	// -- runtype
	189	else if ( !argument.CompareTo("-runtype") ) {
	190	if ((bMissingParam=(++iter>=argc))) break;
	191
	192	TString parameter(argv[iter]);
	193	parameter.Remove(TString::kLeading, ' ');
	194
	195	Bool_t found = kFALSE;
	196
	197	if ( parameter.IsDigit() ) {
	198	Int_t idx = parameter.Atoi();
	199	if ( idx < kRunMax ) {
	200	fRunType = (PprRun_t) idx;
201	found = kTRUE;
202	}
203	}
204	else {
205	for ( Int_t iRun = 0; iRun < kRunMax && !found; iRun++ ) {
206	if ( ! parameter.CompareTo(fgkPprRunName[iRun]) ) {
207	fRunType = (PprRun_t) iRun;
208	found = kTRUE;
209	}
210	}
211	}
212	if ( !found ) {
213	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
214	iResult=-EINVAL;
215	}
216	}
217
218	// -- nevents
219	else if ( !argument.CompareTo("-nevents") ) {
220	if ((bMissingParam=(++iter>=argc))) break;
221
222	TString parameter(argv[iter]);
223	parameter.Remove(TString::kLeading, ' ');
224
225	if ( parameter.IsDigit() ) {
226	fNEvents = parameter.Atoi();
227	}
228
229	else {
230	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
231	iResult=-EINVAL;
232	}
233	}
234
235	// -- coneRadius
236	else if ( !argument.CompareTo("-coneRadius") ) {
237	if ((bMissingParam=(++iter>=argc))) break;
238
239	TString parameter(argv[iter]);
240	parameter.Remove(TString::kLeading, ' ');
241
242	if ( parameter.IsFloat() ) {
243	fJetConeRadius = parameter.Atof();
244	}
245	else {
246	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
247	iResult=-EINVAL;
248	}
249	}
250
251	// -- jetCutMinEt
252	else if ( !argument.CompareTo("-jetCutMinEt") ) {
253	if ((bMissingParam=(++iter>=argc))) break;
254
255	TString parameter(argv[iter]);
256	parameter.Remove(TString::kLeading, ' ');
257
258	if ( parameter.IsFloat() ) {
259	fJetEtMin = parameter.Atof();
260	}
261	else {
262	HLTError("Wrong parameter %s for argument %s.", parameter.Data(), argument.Data());
263	iResult=-EINVAL;
264	}
265	}
266
396e75ed	267	// -- applyParticleCuts
	268	else if ( !argument.CompareTo("-applyParticleCuts") ) {
	269	fApplyParticleCuts = kTRUE;
	270	}
	271
7aac8168	272	// -- Argument not known
	273	else {
	274	HLTError("Unknown argument %s.", argument.Data());
	275	iResult = -EINVAL;
	276	}
	277
	278	} // for ( Int iter = 0; iter<argc && iResult>=0; iter++) {
	279
	280	// -- Check if parameter is missing
	281	if ( bMissingParam ) {
	282	HLTError("Missing parameter for argument %s.", argument.Data());
	283	iResult=-EINVAL;
	284	}
	285
	286	// -- Check if seed and nevents has been specified
	287	if ( !fSeed \|\| !fNEvents ) {
	288	HLTError("-seed and -nevents need to be specified.");
	289	iResult=-EINVAL;
	290	}
	291
	292	// -- Set Random Number seed
	293	gRandom->SetSeed( fSeed );
	294
	295	// -- Run loader
	296	// ---------------
	297	AliPDG::AddParticlesToPdgDataBase();
	298
	299	fRunLoader = AliRunLoader::Open("galice.root","FASTRUN","recreate");
	300	if (fRunLoader == NULL ) {
	301	// gAlice->Fatal("AliHLTMCGeneratorComponent","Can not instatiate the Run Loader");
	302	HLTFatal("Can not instatiate the Run Loader");
	303	return -EINPROGRESS;
	304	}
	305
	306	fRunLoader->SetCompressionLevel(2);
	307	fRunLoader->SetNumberOfEventsPerFile(fNEvents);
	308	fRunLoader->LoadKinematics("RECREATE");
	309	fRunLoader->MakeTree("E");
	310	gAlice->SetRunLoader(fRunLoader);
	311
	312	// -- Create stack
	313	fRunLoader->MakeStack();
	314
	315	// -- Create and Initialize Generator
	316	fGenerator = GeneratorFactory();
	317	fGenerator->SetStack(fRunLoader->Stack());
	318	fGenerator->Init();
	319
	320	HLTInfo("MC Generator setup with : %s", fComment.Data() );
	321
	322	return iResult;
	323	}
	324
	325	// #################################################################################
	326	Int_t AliHLTMCGeneratorComponent::DoDeinit() {
	327	// see header file for class documentation
	328
	329	if ( fpMC )
	330	delete fpMC;
	331	fpMC = NULL;
	332
	333	if ( fpHLTMC )
	334	delete fpHLTMC;
	335	fpHLTMC = NULL;
336
337	return 0;
338	}
339
340	// #################################################################################
341	Int_t AliHLTMCGeneratorComponent::GetEvent( const AliHLTComponentEventData& /evtData/,
342	AliHLTComponentTriggerData& /trigData/,
343	AliHLTUInt8_t* /outputPtr/,
344	AliHLTUInt32_t& size,
345	vector<AliHLTComponentBlockData>& /outputBlocks/ ) {
346	// see header file for class documentation
347
348	if ( !IsDataEvent() )
349	return 0;
350
351	Int_t iResult=0;
352	size=0;
353
354	AliStack* stack = fRunLoader->Stack();
355	AliHeader* header = fRunLoader->GetHeader();
356
357	// -- Initialize event
358	header->Reset(0,fCurrentEvent);
359
360	// -- Reset stack
361	stack->Reset();
362
363	// -- Generate event
364	fGenerator->Generate();
365
366	// -- Finish event
367	header->SetNprimary(stack->GetNprimary());
368	header->SetNtrack(stack->GetNtrack());
369
370	// -- Finish Event
371	stack->FinishEvent();
372	header->SetStack(stack);
373
374	// -- Create HLT MC Event
375	if ( fpHLTMC )
376	delete fpHLTMC;
396e75ed	377	fpHLTMC = new AliHLTMCEvent( fApplyParticleCuts);
7aac8168	378
	379	// -- Fill HLT MC event
	380	iResult = fpHLTMC->FillMCEvent( stack, header );
	381
	382	// -- Send HLT MC event
	383	if ( fpHLTMC && !iResult )
	384	PushBack( fpHLTMC, kAliHLTDataTypeMCObject\|kAliHLTDataOriginHLT, 0 );
	385
	386	fCurrentEvent++;
	387
	388	return iResult;
	389	}
	390
	391	// #################################################################################
	392	AliGenerator* AliHLTMCGeneratorComponent::GeneratorFactory() {
	393
	394	// -- Configuration for hard QCD processes generation with PYTHIA
	395	// ----------------------------------------------------------------
	396	AliGenPythia * gener = new AliGenPythia(-1);
	397
	398	gener->SetEnergyCMS(fEcms); // Centre of mass energy
	399	gener->SetGluonRadiation(1,1);
	400	gener->SetForceDecay(kAll); // Decay type
	401
	402	// -- Final state kinematic cuts
	403	gener->SetJetEtaRange(fJetEtaMin, fJetEtaMax);
	404	gener->SetJetPhiRange(0., 360.);
	405	gener->SetJetEtRange(fJetEtMin, fJetEtMax);
	406	// gener->SetMomentumRange(0,99999999);
	407
	408	// -- Structure function
	409	gener->SetStrucFunc(kCTEQ4L);
	410
	411	/************************************************************************
	412	* void AliGenPythia::SetPycellParameters()
	413	* Float_t etamax -- fPycellEtaMax -- PARU(51) -- +/- eta max of detector
	414	* Int_t neta -- fPycellNEta -- MSTU(51) -- N cells in eta
	415	* Int_t nphi -- fPycellNPhi -- MSTU(52) -- N cells in phi
	416	* Float_t thresh -- fPycellThreshold -- PARU(58) -- cells with Et below treshold are discared
	417	* Float_t etseed -- fPycellEtSeed -- PARU(52) -- Et threshold for seed cells
	418	* Float_t minet -- fPycellMinEtJet -- PARU(53) -- Min Et for jets
	419	* Float_t r -- fPycellMaxRadius -- PARU(54) -- cone radius
	420	************************************************************************/
	421	gener->SetPycellParameters(fJetEtaMax, 274, 432, 0., 4., fJetEtMin, fJetConeRadius);
	422
	423	// gener->SetPtKick(5.); // set the intrinsic kt to 5 GeV/c
	424
	425	// -- Jet Quenching
	426	gener->SetQuench(fQuenching);
	427	if( fQuenching == 1){
	428	Float_t k = 6e5*(fQhat/1.7) ; // qhat=1.7, k = 6e5, default value
	429	AliPythia::Instance()->InitQuenching(0.,0.1,k,0,0.95,6);
	430	}
	431
	432	switch (fRunType) {
	433	case kPythia6Jets20_24:
	434	{
	435	fComment = fComment.Append(":Pythia jets 20-24 GeV @ 14 TeV");
	436	gener->SetProcess(kPyJets);// Process type
	437	gener->SetPtHard(20., 24.);// Pt transfer of the hard scattering
	438	fPtHardMin=20.;
	439	}
	440	break;
	441	case kPythia6Jets24_29:
442	{
443	fComment = fComment.Append(":Pythia jets 24-29 GeV @ 14 TeV");
444	gener->SetProcess(kPyJets);// Process type
445	gener->SetPtHard(24., 29.);// Pt transfer of the hard scattering
446	fPtHardMin=24.;
447	}
448	break;
449	case kPythia6Jets29_35:
450	{
451	fComment = fComment.Append(":Pythia jets 29-35 GeV @ 14 TeV");
452	gener->SetProcess(kPyJets);// Process type
453	gener->SetPtHard(29., 35.);// Pt transfer of the hard scattering
454	fPtHardMin=29.;
455	}
456	break;
457	case kPythia6Jets35_42:
458	{
459	fComment = fComment.Append(":Pythia jets 35-42 GeV @ 14 TeV");
460	gener->SetProcess(kPyJets);// Process type
461	gener->SetPtHard(35., 42.);// Pt transfer of the hard scattering
462	fPtHardMin=35.;
463	}
464	break;
465	case kPythia6Jets42_50:
466	{
467	fComment = fComment.Append(":Pythia jets 42-50 GeV @ 14 TeV");
468	gener->SetProcess(kPyJets);// Process type
469	gener->SetPtHard(42., 50.);// Pt transfer of the hard scattering
470	fPtHardMin=42.;
471	}
472	break;
473	case kPythia6Jets50_60:
474	{
475	fComment = fComment.Append(":Pythia jets 50-60 GeV @ 14 TeV");
476	gener->SetProcess(kPyJets);// Process type
477	gener->SetPtHard(50., 60.);// Pt transfer of the hard scattering
478	fPtHardMin=50.;
479	}
480	break;
481	case kPythia6Jets60_72:
482	{
483	fComment = fComment.Append(":Pythia jets 60-72 GeV @ 14 TeV");
484	gener->SetProcess(kPyJets);// Process type
485	gener->SetPtHard(60., 72.);// Pt transfer of the hard scattering
486	fPtHardMin=60.;
487	}
488	break;
489	case kPythia6Jets72_86:
490	{
491	fComment = fComment.Append(":Pythia jets 72-86 GeV @ 14 TeV");
492	gener->SetProcess(kPyJets);// Process type
493	gener->SetPtHard(72., 86.);// Pt transfer of the hard scattering
494	fPtHardMin=72.;
495	}
496	break;
497	case kPythia6Jets86_104:
498	{
499	fComment = fComment.Append(":Pythia jets 86-104 GeV @ 14 TeV");
500	gener->SetProcess(kPyJets);// Process type
501	gener->SetPtHard(86., 104.);// Pt transfer of the hard scattering
502	fPtHardMin=86.;
503	}
504	break;
505	case kPythia6Jets104_125:
506	{
507	fComment = fComment.Append(":Pythia jets 104-125 GeV @ 14 TeV");
508	gener->SetProcess(kPyJets);// Process type
509	gener->SetPtHard(104., 125.);// Pt transfer of the hard scattering
510	fPtHardMin=104.;
511	}
512	break;
513	case kPythia6Jets125_150:
514	{
515	fComment = fComment.Append(":Pythia jets 125-150 GeV @ 14 TeV");
516	gener->SetProcess(kPyJets);// Process type
517	gener->SetPtHard(125., 150.);// Pt transfer of the hard scattering
518	fPtHardMin=125.;
519	}
520	break;
521	case kPythia6Jets150_180:
522	{
523	fComment = fComment.Append(":Pythia jets 150-180 GeV @ 14 TeV");
524	gener->SetProcess(kPyJets);// Process type
525	gener->SetPtHard(150., 180.);// Pt transfer of the hard scattering
526	fPtHardMin=150.;
527	}
528	break;
529	case kPyJetJet:
530	{
531	fComment = fComment.Append(" Jet-jet @ 14 TeV");
532	gener->SetProcess(kPyJets);
533	gener->SetPtHard(fPtHardMin,fPtHardMax);
534	gener->SetEventListRange(0,1); // XXXX
535	}
536	break;
537	case kPyGammaJetPHOS:
538	{
539	fComment = fComment.Append(" pp->jet + Gamma over PHOS @ 14 TeV");
540	gener->SetProcess(kPyDirectGamma);
541	gener->SetPtHard(fPtHardMin,fPtHardMax);
542	//gener->SetYHard(-1.0,1.0);
543	gener->SetGammaEtaRange(-0.13,0.13);
544	gener->SetGammaPhiRange(218.,322.); //Over 5 modules +-2 deg
545	gener->SetEventListRange(0,1); // XXXX
546	}
547	break;
548
549	default: break;
550	}
551
552	return gener;
553	}