PWGJE/macros/examples/pythia6FastjetExample.C

   1 void run(Int_t nEvent = 50, Float_t e_cms = 2760) {
   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("libSISConePlugin");
  13   gSystem->Load("libJETANdev");
  14   gSystem->Load("libFASTJETANdev");
  15
  16   gSystem->Load("libpythia6.so");
  17   gSystem->Load("libEGPythia6.so");
  18   gSystem->Load("libAliPythia6.so");
  19
  20
  21   AliPDG::AddParticlesToPdgDataBase(); // to add some PDF codes to TDatabasePDG
  22
  23   // Create random number generator and set seed
  24   AliPythiaRndm::SetPythiaRandom(new TRandom3());
  25   AliPythiaRndm::GetPythiaRandom()->SetSeed(clock()+gSystem->GetPid());
  26
  27   AliPythia6 *pythia=AliPythia6::Instance();
  28
  29   pythia->SetCKIN(3,10);   // minimum hard pt
  30   pythia->SetCKIN(4,1000);  // maximum hard pt
  31
  32   pythia->SetMDCY(pythia->Pycomp(111),1,0);  // switch off pi0 decay
  33
  34   pythia->Initialize("CMS","p","p",e_cms);
  35
  36
  37   AliFastJetHeaderV1 *header = new AliFastJetHeaderV1;
  38   header->SetBGMode(0);
  39   header->SetRadius(0.4);
  40   //header->SetGhostEtaMax(2);
  41   //header->SetGhostArea(0.05);
  42   header->SetAlgorithm(2); // antikt_algorithm = 2, kt = 0 (see fastjet/fastjet/JetDefinition.hh
  43
  44   AliFastJetFinder *FastJet = new AliFastJetFinder;
  45   FastJet->SetJetHeader(header);
  46
  47   AliAODEvent *aod = new AliAODEvent();
  48   aod->CreateStdContent();
  49   FastJet->ConnectAOD(aod);
  50
  51   AliJetCalTrkEvent JetFinderEvent(0,1);
  52   TClonesArray *plist = new TClonesArray("TParticle");
  53   TClonesArray aliplist("AliMCParticle",1000);
  54
  55   for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
  56
  57     TProcessID::SetObjectCount(0); // Needed for TRefs in AliCalTrkTrack and AliAODJet
  58
  59     pythia->GenerateEvent();
  60
  61     pythia->GetParticles(plist);
  62
  63     aliplist.Clear();
  64     JetFinderEvent.Clear();
  65
  66     Int_t n_part = plist->GetEntries();
  67     for (Int_t i_part = 0; i_part < n_part; i_part++) {
  68       part=(TParticle*) plist->At(i_part);
  69
  70       if (part->GetStatusCode() >= 10)  // Not a final state particle
  71         continue;
  72
  73       new (aliplist[i_part]) AliMCParticle(part);
  74
  75       JetFinderEvent.AddCalTrkTrackKine((AliMCParticle*)aliplist[i_part],1,1);
  76     }
  77
  78     aod->ClearStd();
  79     FastJet->Reset();
  80     FastJet->SetCalTrkEvent(JetFinderEvent);
  81     FastJet->ProcessEvent();
  82     if (aod->GetNJets() > 0) {
  83       cout << "event " << iEvent << " " << aod->GetNJets() << " jets found" << endl;
  84       for (Int_t iJet = 0; iJet < aod->GetNJets(); iJet++) {
  85         jet = aod->GetJet(iJet);
  86         cout << "\t jet " << iJet << " pt " << jet->Pt() << " eta " << jet->Eta() << " phi " << jet->Phi() << endl;
  87       }
  88     }
  89   }
  90 }