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("libJETAN");
  14   gSystem->Load("libFASTJETAN");
  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->SetRparam(0.4);
  41   //header->SetGhostEtaMax(2);
  42   //header->SetGhostArea(0.05);
  43   header->SetAlgorithm(2); // antikt_algorithm = 2, kt = 0 (see fastjet/fastjet/JetDefinition.hh
  44
  45   AliFastJetFinder *FastJet = new AliFastJetFinder;
  46   FastJet->SetJetHeader(header);
  47
  48   AliAODEvent *aod = new AliAODEvent();
  49   aod->CreateStdContent();
  50   FastJet->ConnectAOD(aod);
  51
  52   AliJetCalTrkEvent JetFinderEvent(0,1);
  53   TClonesArray *plist = new TClonesArray("TParticle");
  54   TClonesArray aliplist("AliMCParticle",1000);
  55
  56   for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
  57
  58     TProcessID::SetObjectCount(0); // Needed for TRefs in AliCalTrkTrack and AliAODJet
  59
  60     pythia->GenerateEvent();
  61
  62     pythia->GetParticles(plist);
  63
  64     aliplist.Clear();
  65     JetFinderEvent.Clear();
  66
  67     Int_t n_part = plist->GetEntries();
  68     for (Int_t i_part = 0; i_part < n_part; i_part++) {
  69       part=(TParticle*) plist->At(i_part);
  70
  71       if (part->GetStatusCode() >= 10)  // Not a final state particle
  72         continue;
  73
  74       new (aliplist[i_part]) AliMCParticle(part);
  75
  76       JetFinderEvent.AddCalTrkTrackKine((AliMCParticle*)aliplist[i_part],1,1);
  77     }
  78
  79     aod->ClearStd();
  80     FastJet->Reset();
  81     FastJet->SetCalTrkEvent(JetFinderEvent);
  82     FastJet->ProcessEvent();
  83     if (aod->GetNJets() > 0) {
  84       cout << "event " << iEvent << " " << aod->GetNJets() << " jets found" << endl;
  85       for (Int_t iJet = 0; iJet < aod->GetNJets(); iJet++) {
  86         jet = aod->GetJet(iJet);
  87         cout << "\t jet " << iJet << " pt " << jet->Pt() << " eta " << jet->Eta() << " phi " << jet->Phi() << endl;
  88       }
  89     }
  90   }
  91 }