PWGJE/macros/examples/pythia6FastjetExample.C

   1 void pythia6FastjetExample(Int_t nEvent = 50, Float_t e_cms = 2760, Float_t pthardmin = 10, Float_t pthardmax = 50) {
   2
   3   // example macro for on-the-fly generation of PYTHIA6 events
   4   // and analysis with new reader interface and FastJet
   5   // M. van Leeuwen
   6
   7   gSystem->Load("libANALYSIS");
   8   gSystem->Load("libANALYSISalice");
   9   gSystem->Load("libCGAL");
  10   gSystem->Load("libfastjet");
  11   gSystem->Load("libsiscone");
  12   gSystem->Load("libsiscone_spherical");
  13   gSystem->Load("libfastjetplugins");
  14   gSystem->Load("libJETAN");
  15   gSystem->Load("libFASTJETAN");
  16
  17   gSystem->Load("libpythia6");
  18   gSystem->Load("libEGPythia6");
  19   gSystem->Load("liblhapdf");
  20   gSystem->Load("libAliPythia6");
  21
  22
  23   AliPDG::AddParticlesToPdgDataBase(); // to add some PDF codes to TDatabasePDG
  24
  25   // Create random number generator and set seed
  26   AliPythiaRndm::SetPythiaRandom(new TRandom3());
  27   AliPythiaRndm::GetPythiaRandom()->SetSeed(clock()+gSystem->GetPid());
  28
  29   AliGenPythia *pythia=new AliGenPythia(1);
  30   pythia->SetProcess(kPyJets);
  31   pythia->SetPtHard(pthardmin, pthardmax);
  32   pythia->SetEnergyCMS(e_cms);
  33
  34   pythia->Init();
  35
  36   //
  37   // Jet finder settings go via the FastJetHeader
  38   //
  39   AliFastJetHeaderV1 *header = new AliFastJetHeaderV1;
  40   header->SetBGMode(0);
  41   //  header->SetRadius(0.4);
  42   header->SetRparam(0.4);
  43   //header->SetGhostEtaMax(2);
  44   //header->SetGhostArea(0.05);
  45   header->SetAlgorithm(2); // antikt_algorithm = 2, kt = 0 (see fastjet/fastjet/JetDefinition.hh
  46
  47   AliFastJetFinder *FastJet = new AliFastJetFinder;
  48   FastJet->SetJetHeader(header);
  49
  50   // Infrastructure needed for AliGenPythia
  51   AliStack *stack = new AliStack();
  52   // There is a big mess with the connection between runloader, AliRun and the gAlice pointer.
  53   // This order of things seems to work...
  54   AliRunLoader *dummyrl = new AliRunLoader("dummyEvtFolder");
  55   dummyrl->MakeHeader();
  56   dummyrl->SetEventNumber(0);
  57   gAlice->SetRunLoader(dummyrl);
  58   pythia->SetStack(stack);
  59
  60   // Set up dummy AOD event for JETAN output
  61   AliAODEvent *aod = new AliAODEvent();
  62   aod->CreateStdContent();
  63   FastJet->ConnectAOD(aod);
  64
  65   // Set up input structures for FastJet
  66   AliJetCalTrkEvent JetFinderEvent(0,1);
  67   TClonesArray aliplist("AliMCParticle",1000);
  68
  69   for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
  70
  71     TProcessID::SetObjectCount(0); // Needed for TRefs in AliCalTrkTrack and AliAODJet
  72
  73     stack->Reset();
  74     pythia->Generate();
  75
  76     aliplist.Clear();
  77     JetFinderEvent.Clear();
  78
  79     Int_t n_part = stack->GetNtrack();
  80     for (Int_t i_part = 0; i_part < n_part; i_part++) {
  81       part=(TParticle*) stack->Particle(i_part);
  82
  83       if (part->GetStatusCode() >= 10)  // Not a final state particle
  84         continue;
  85
  86       new (aliplist[i_part]) AliMCParticle(part);
  87
  88       JetFinderEvent.AddCalTrkTrackKine((AliMCParticle*)aliplist[i_part],1,1);
  89     }
  90
  91     aod->ClearStd();
  92     FastJet->Reset();
  93     FastJet->SetCalTrkEvent(JetFinderEvent);
  94     FastJet->ProcessEvent();
  95     if (aod->GetNJets() > 0) {
  96       cout << "event " << iEvent << " " << aod->GetNJets() << " jets found" << endl;
  97       for (Int_t iJet = 0; iJet < aod->GetNJets(); iJet++) {
  98         jet = aod->GetJet(iJet);
  99         cout << "\t jet " << iJet << " pt " << jet->Pt() << " eta " << jet->Eta() << " phi " << jet->Phi() << endl;
 100       }
 101     }
 102   }
 103 }
 104
 105 // kept for backward compatibility; this was the old function name
 106 void run(Int_t nEvent = 50, Float_t e_cms = 2760) {
 107   pythia6FastjetExample(nEvent, e_cms);
 108 }