PWGUD/macros/genLevelSimulation/fastGen.C

   1 typedef enum {kPhojet = -1,kPyTuneCDFA=100,kPyTuneAtlasCSC=306, kPyTuneCMS6D6T=109, kPyTunePerugia0=320 } Tune_t;
   2
   3 AliGenerator*  CreateGenerator(Tune_t tune = kPyTuneCDFA , Float_t energy);
   4
   5 void fastGen(Tune_t tune = kPyTuneCDFA , Float_t energy, Int_t nev = 1, TString process)
   6 {
   7   // Add all particles to the PDG database
   8   AliPDG::AddParticlesToPdgDataBase();
   9
  10   // set the random seed
  11   TDatime date;
  12   UInt_t seed    = date.Get()+gSystem->GetPid();
  13   gRandom->SetSeed(seed);
  14   cout<<"Seed for random number generation= "<<seed<<endl;
  15
  16
  17   //  Runloader
  18   AliRunLoader* rl = AliRunLoader::Open("galice.root", "FASTRUN","recreate");
  19
  20   rl->SetCompressionLevel(2);
  21   rl->SetNumberOfEventsPerFile(nev);
  22   rl->LoadKinematics("RECREATE");
  23   rl->MakeTree("E");
  24   gAlice->SetRunLoader(rl);
  25
  26   //  Create stack
  27   rl->MakeStack();
  28   AliStack* stack      = rl->Stack();
  29
  30   //  Header
  31   AliHeader* header = rl->GetHeader();
  32   //
  33   //  Create and Initialize Generator
  34   AliGenerator *gener = CreateGenerator(tune,energy);
  35   gener->Init();
  36   // if nsd switch off single diffraction
  37   if ( process == "NSD"){
  38     if(tune != kPhojet) {
  39       AliPythia::Instance()->   SetMSUB(92,0);             // single diffraction AB-->XB
  40       AliPythia::Instance()-> SetMSUB(93,0);             // single diffraction AB-->AX
  41     }
  42     else {
  43       cout << "NSD not yet implemented in the phojet case" << endl;
  44       exit(1);
  45     }
  46   }
  47   gener->SetStack(stack);
  48
  49   //
  50   //                        Event Loop
  51   //
  52   Int_t iev;
  53
  54   for (iev = 0; iev < nev; iev++) {
  55
  56     if(!(iev%500)) printf("\n \n Event number %d \n \n", iev);
  57
  58     //  Initialize event
  59     header->Reset(0,iev);
  60     rl->SetEventNumber(iev);
  61     stack->Reset();
  62     rl->MakeTree("K");
  63     //  stack->ConnectTree();
  64
  65     //  Generate event
  66     gener->Generate();
  67     //  Analysis
  68     //  Int_t npart = stack->GetNprimary();
  69     //  printf("Analyse %d Particles\n", npart);
  70     //  for (Int_t part=0; part<npart; part++) {
  71     //      TParticle *MPart = stack->Particle(part);
  72     //      Int_t mpart  = MPart->GetPdgCode();
  73     //      printf("Particle %d\n", mpart);
  74     //  }
  75
  76     //  Finish event
  77     header->SetNprimary(stack->GetNprimary());
  78     header->SetNtrack(stack->GetNtrack());
  79     //      I/O
  80     //
  81     stack->FinishEvent();
  82     header->SetStack(stack);
  83     rl->TreeE()->Fill();
  84     rl->WriteKinematics("OVERWRITE");
  85
  86   } // event loop
  87     //
  88     //                         Termination
  89     //  Generator
  90   gener->FinishRun();
  91   //  Write file
  92   rl->WriteHeader("OVERWRITE");
  93   gener->Write();
  94   rl->Write();
  95
  96 }
  97
  98
  99 AliGenerator*  CreateGenerator(Tune_t tune, Float_t energy)
 100 {
 101
 102
 103
 104   if (tune == -1) {
 105
 106     // phojet
 107     AliGenDPMjet* gener = new AliGenDPMjet(1);
 108     gener->SetProcess(kDpmMb);
 109     gener->SetProjectile("P", 1, 1);
 110     gener->SetTarget("P", 1, 1);
 111
 112     gener->SetEnergyCMS(energy);
 113     return gener;
 114   }
 115
 116   if (tune != kPyTuneAtlasCSC && tune != kPyTuneCDFA && tune != kPyTuneCMS6D6T && tune != kPyTunePerugia0) {
 117
 118     Printf("Unknown pythia tune, quitting");
 119     exit(1);
 120
 121   } else {
 122
 123     AliGenPythia * gener =  new AliGenPythia(1);
 124     //
 125     //  set pythia tune
 126     gener->SetTune(tune);
 127
 128     //   structure function
 129     if(tune == kPyTuneAtlasCSC) {
 130       cout << "Setting structure function" << endl;
 131       gener->SetStrucFunc(kCTEQ61);
 132     }
 133     if(tune == kPyTuneCMS6D6T) {
 134       cout << "Setting structure function" << endl;
 135       gener->SetStrucFunc(kCTEQ6l);
 136     }
 137     if(tune == kPyTunePerugia0) {
 138       cout << "Setting new parton shower" << endl;
 139       gener->UseNewMultipleInteractionsScenario();
 140     }
 141     //   charm, beauty, charm_unforced, beauty_unforced, jpsi, jpsi_chi, mb
 142     gener->SetProcess(kPyMb);
 143     //   Centre of mass energy
 144     gener->SetEnergyCMS(energy);
 145
 146     // Set target/projectile // Is this needded?
 147     gener->SetProjectile("P", 1, 1);
 148     gener->SetTarget("P", 1, 1);
 149
 150     return gener;
 151   }
 152 }
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164