1 void run(Int_t nEvent = 50, Float_t e_cms = 2760) {
3 // example macro for on-the-fly generation of PYTHIA6 events
4 // and analysis with new reader interface and FastJet
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");
16 gSystem->Load("libpythia6.so");
17 gSystem->Load("libEGPythia6.so");
18 gSystem->Load("libAliPythia6.so");
21 AliPDG::AddParticlesToPdgDataBase(); // to add some PDF codes to TDatabasePDG
23 // Create random number generator and set seed
24 AliPythiaRndm::SetPythiaRandom(new TRandom3());
25 AliPythiaRndm::GetPythiaRandom()->SetSeed(clock()+gSystem->GetPid());
27 AliPythia6 *pythia=AliPythia6::Instance();
29 pythia->SetCKIN(3,10); // minimum hard pt
30 pythia->SetCKIN(4,1000); // maximum hard pt
32 pythia->SetMDCY(pythia->Pycomp(111),1,0); // switch off pi0 decay
34 pythia->Initialize("CMS","p","p",e_cms);
37 AliFastJetHeaderV1 *header = new AliFastJetHeaderV1;
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
45 AliFastJetFinder *FastJet = new AliFastJetFinder;
46 FastJet->SetJetHeader(header);
48 AliAODEvent *aod = new AliAODEvent();
49 aod->CreateStdContent();
50 FastJet->ConnectAOD(aod);
52 AliJetCalTrkEvent JetFinderEvent(0,1);
53 TClonesArray *plist = new TClonesArray("TParticle");
54 TClonesArray aliplist("AliMCParticle",1000);
56 for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
58 TProcessID::SetObjectCount(0); // Needed for TRefs in AliCalTrkTrack and AliAODJet
60 pythia->GenerateEvent();
62 pythia->GetParticles(plist);
65 JetFinderEvent.Clear();
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);
71 if (part->GetStatusCode() >= 10) // Not a final state particle
74 new (aliplist[i_part]) AliMCParticle(part);
76 JetFinderEvent.AddCalTrkTrackKine((AliMCParticle*)aliplist[i_part],1,1);
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;