HLT/exa/MakePythia.C

   1 //$Id$
   2
   3 /* A little macro filling a root
   4    tree with PYTHIA6 events. Same
   5    as generate_PYTHIA.C */
   6
   7 typedef struct
   8 {
   9   Int_t    no;
  10   Int_t    part;
  11   Int_t    total;
  12   Int_t    npart;
  13   Int_t    nenergy;
  14   Int_t    nhard;
  15   Float_t  b;
  16   Float_t  phi;
  17 } THeader;
  18
  19 UInt_t makeSeed(int mode=0)
  20 {
  21   switch (mode) {
  22   case 1:
  23     {
  24       TDatime date;
  25       return date.GetDate();
  26     }
  27   case 2:
  28     {
  29       TDatime date;
  30       UInt_t  seed1 = (date.GetDate() / 1000000 * 100 +
  31                        date.GetTime() / 1000000);
  32       UInt_t  seed2 = (date.GetTime() % 10000);
  33       TRanMar ranmar(seed1, seed2);
  34       Int_t   eat = (date.GetDate() - 19980101 + seed2 +
  35                      gSystem->GetPid());
  36       for (Int_t i = 0; i < (eat + 250000) ; i++)
  37         ranmar.Rndm();
  38       return 2 * Int_t(ranmar.Rndm() * (TMath::Power(2,30) - 1)) + 1;
  39     }
  40   case 3:
  41     {
  42       TDatime date;
  43       UInt_t  seed1 = (date.GetDate() / 1000000 * 100 +
  44                        date.GetTime() / 1000000 + gSystem->GetPid());
  45       TRandom rand(seed1);
  46       return Int_t(rand.Rndm() * (TMath::Power(2,30) - 1)) + 1;
  47     }
  48
  49     break;
  50   }
  51   return 0;
  52 }
  53
  54 Int_t MakePythia(Int_t nEvents=10,Char_t *rfile="pythia6.root")
  55 {
  56    gSystem->Load("$ROOTSYS/lib/libPythia6.so");
  57    gSystem->Load("$ROOTSYS/lib/libEG.so");                   // Root Event Generator interface
  58    gSystem->Load("$ROOTSYS/lib/libEGPythia6.so");            // Root interface to Pythia6
  59    //gDebug=1;
  60
  61    TStopwatch*   stopWatch = new TStopwatch();
  62    TPythia6* pythia6       = new TPythia6();
  63    TClonesArray* particles = new TClonesArray("TParticle");
  64    THeader header;
  65
  66    Char_t filename[1024];
  67    sprintf(filename,"%s",rfile);
  68
  69    TFile* file = new TFile(filename,"RECREATE");
  70    file->SetCompressionLevel(1);
  71    TTree* tree = new TTree("pythia","pythia events");
  72    tree->Branch("particles",&particles);
  73    //tree->Branch("header",&header, "no/I:part:total:npart:nenergy:nhard:b/F:phi");
  74
  75    stopWatch->Start();
  76
  77    //      select Pythia min. bias model, taken from AliPythia.C
  78    pythia6->SetMSEL(0);
  79    pythia6->SetMSUB(92,1);      // single diffraction AB-->XB
  80    pythia6->SetMSUB(93,1);      // single diffraction AB-->AX
  81    pythia6->SetMSUB(94,1);      // double diffraction
  82    pythia6->SetMSUB(95,1);      // low pt production
  83    pythia6->SetMSTP(81,1);      // multiple interactions switched on
  84    pythia6->SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
  85    pythia6->SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of
  86    // multiple interaction at a reference energy = 14000 GeV
  87    pythia6->SetPARP(89,14000.); // reference energy for the above parameter
  88    pythia6->SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
  89
  90    // set fragmentation on (default)
  91    pythia6->SetMSTP(111,1);
  92
  93    // don't smear the primary vertex
  94    pythia6->SetMSTP(151,0);
  95
  96    //pythia6->SetMSTP(61,0);
  97    //pythia6->SetMSTP(64,0);
  98
  99    pythia6->SetMRPY(1,makeSeed(3));
 100    pythia6->Initialize("cms","p","p",5500);
 101    //pythia6->Initialize("cms","p","p",1800);
 102
 103    stopWatch->Stop();
 104    stopWatch->Print();
 105
 106    Int_t i;
 107    for (i = 0; i < nEvents; i++) {
 108      stopWatch->Start();
 109      cout << "Event # " << i << " ... " << flush;
 110
 111      pythia6->GenerateEvent();
 112
 113      header.no    = i;
 114      header.total = pythia6->ImportParticles(particles,"All");
 115      header.part  = 0;
 116
 117      header.npart    = 0;
 118      header.nenergy  = 0;
 119      header.nhard    = 0;
 120
 121      header.b     = 0;
 122      header.phi   = 0;
 123
 124      //particles->Print();
 125      //pythia6->Pylist(1);
 126
 127      tree->Fill();
 128
 129      cout << "done (" << header.total << " particles) " << flush;
 130      stopWatch->Stop();
 131      stopWatch->Print();
 132    }
 133
 134    file->Write();
 135    file->Close();
 136
 137   return 0;
 138 }
 139