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   // 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 }