MUON/fastMUONGen.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /// \ingroup macros
  19 /// \file fastMUONGen.C
  20 /// \brief An example how to use the AliGenMUONCocktailpp generator
  21 /// without GEANT simulation of the detector
  22 ///
  23 /// \author H. Woehri,  A. De Falco, INFN Cagliari, April 2007
  24 ///
  25 /// The macro switches on all decay modes for the resonances,
  26 /// while for the minimum bias Pythia event we additionally
  27 /// switch on the muonic decays of pions and Kaons
  28 /// Its outcome can be further processed by the macro
  29 /// fastMUONSim.C
  30
  31 #if !defined(__CINT__) || defined(__MAKECINT__)
  32 #include "AliGenerator.h"
  33 #include "AliRunLoader.h"
  34 #include "AliRun.h"
  35 #include "AliHeader.h"
  36 #include "AliStack.h"
  37 #include "AliGenCocktail.h"
  38 #include "AliDecayerPythia.h"
  39 #include "AliPDG.h"
  40 #include "AliGenMUONCocktailpp.h"
  41
  42 #include "TParticle.h"
  43 #endif
  44
  45 AliGenerator*  CreateGeneratorMC(Int_t mult);
  46 Int_t SetupOutputDirectory();
  47
  48 /// The third argument \em mult allows to select at a "pre-trigger" level
  49 /// events with "mult" muons in the MUON detector's phase space
  50 /// window, which is defined in the method "CreateGeneratorMC(Int_t mult)".
  51 /// Note that in this routine also a cut on the muon's origin
  52 /// should be placed in order not to trigger on muons from pi/K decays
  53 /// that were decayed by Pythia way inside the absorber or in the muon
  54 /// spectrometer itself
  55 void fastMUONGen(Int_t nev = 1, char* filename = "galice.root", Int_t mult = 2)
  56 {
  57   Int_t runNumber = SetupOutputDirectory();
  58   gAlice->SetRunNumber(runNumber);
  59   printf("\n\n\n\nsetting run number to %d\n\n\n\n", runNumber);
  60
  61   //  Update data base
  62   AliPDG::AddParticlesToPdgDataBase();
  63   //
  64   AliRunLoader* rl = AliRunLoader::Open(filename,"FASTRUN","recreate");
  65   rl->SetCompressionLevel(2);
  66   rl->SetNumberOfEventsPerFile(10000);
  67   rl->LoadKinematics("RECREATE");
  68   rl->MakeTree("E");
  69   gAlice->SetRunLoader(rl);
  70   //
  71   rl->MakeStack();
  72   AliStack* stack      = rl->Stack();
  73   AliHeader* header = rl->GetHeader();
  74   //
  75   //  Create and Initialize Generator
  76   AliGenerator *gener = CreateGeneratorMC(mult);
  77   gener->SetStack(stack);
  78   gener->Init();
  79   //
  80   for (int iev = 0; iev < nev; iev++) {
  81
  82     if(iev%1000 == 0) printf("Event %d\n", iev);
  83     //  Initialize event
  84     header->Reset(0,iev);
  85     rl->SetEventNumber(iev);
  86     stack->Reset();
  87     rl->MakeTree("K");
  88
  89     gener->Generate();
  90 /*
  91     // Printing
  92     Int_t npart = stack->GetNprimary();
  93     for (Int_t part=0; part<npart; part++) {
  94       TParticle *MPart = stack->Particle(part);
  95       Int_t mpart  = MPart->GetPdgCode();
  96       printf("Particle %5d %5d  %5d\n", part, mpart, MPart->GetFirstMother());
  97     }
  98 */
  99
 100     //  Finish event
 101     header->SetNprimary(stack->GetNprimary());
 102     header->SetNtrack(stack->GetNtrack());
 103     //      I/O
 104     //
 105     stack->FinishEvent();
 106     header->SetStack(stack);
 107     rl->TreeE()->Fill();
 108     rl->WriteKinematics("OVERWRITE");
 109   } // event loop
 110
 111   gener->FinishRun();
 112   rl->WriteHeader("OVERWRITE");
 113   gener->Write();
 114   rl->Write();
 115 }
 116 //===============================================
 117 AliGenerator*  CreateGeneratorMC(Int_t mult)
 118 {
 119     AliGenMUONCocktailpp* gener = new AliGenMUONCocktailpp();
 120     gener->SetPtRange(0.,100.);
 121     gener->SetYRange(-4.,-2.4);
 122     gener->SetPhiRange(0., 360.);
 123     gener->SetMuonMultiplicity(mult);
 124     gener->SetMuonPtCut(0.5);
 125     gener->SetMuonThetaRange(171.,178.);
 126     gener->SetOrigin(0.,0.,0.);
 127     gener->SetSigma(0.,0.,5.);
 128     gener->SetVertexSmear(kPerEvent);
 129     //for resonances: "kAll" all decay modes are openend
 130     gener->SetDecayModeResonance(kAll);
 131     //for MC Pythia events: all decays, including muonic decays from pion and Kaons are opened
 132     gener->SetDecayModePythia(kAllMuonic);
 133     gener->SetMuonOriginCut(-130.); //prevent the trigger muon(s) from decaying after 130 cm
 134     AliDecayerPythia* decayer = new AliDecayerPythia();
 135     gener->SetDecayer(decayer);
 136     return gener;
 137 }
 138 //===============================================
 139 Int_t SetupOutputDirectory(){
 140
 141   //Setting up the name for the output directory:
 142   static Int_t sseed = 0; // Set 0 to use the current time
 143   gRandom->SetSeed(sseed);
 144   UInt_t theSeed = gRandom->GetSeed();
 145   Int_t labIndex = 5; //1 subatech, 2 clermont, 3 Torino, 4 Orsay
 146                       // 5 Cagliari, etc...
 147   Int_t runNumber = (theSeed%100000000 + 100000000*labIndex);
 148
 149   Char_t name[100];
 150   sprintf(name, "touch run_%d", runNumber);
 151   gSystem->Exec(name);
 152
 153   return runNumber;
 154 }