[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.) {
  // 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( ULong_t& constBase, Double_t& 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;
      }
    } 

    // -- 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();

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