HLT/BASE/util/AliHLTMCGeneratorComponent.cxx

   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   fApplyParticleCuts(kFALSE) {
  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 // #################################################################################
 136 void AliHLTMCGeneratorComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) {
 137   // see header file for class documentation
 138
 139   constBase = 50000;
 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
 267     // -- applyParticleCuts
 268     else if ( !argument.CompareTo("-applyParticleCuts") ) {
 269       fApplyParticleCuts = kTRUE;
 270     }
 271
 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;
 377   fpHLTMC = new AliHLTMCEvent( fApplyParticleCuts);
 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 }