]>
Commit | Line | Data |
---|---|---|
c45ff98f | 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 | AliPythia6 *pythia=AliPythia6::Instance(); | |
24 | ||
25 | pythia->SetCKIN(3,10); // minimum hard pt | |
26 | pythia->SetCKIN(4,1000); // maximum hard pt | |
27 | ||
28 | pythia->SetMRPY(1,clock()+gSystem->GetPid()); // set random seed | |
29 | pythia->SetMRPY(2,0); // Set to 0 to pick up seed | |
30 | ||
31 | pythia->SetMDCY(pythia->Pycomp(111),1,0); // switch off pi0 decay | |
32 | ||
33 | pythia->Initialize("CMS","p","p",e_cms); | |
34 | ||
35 | ||
36 | AliFastJetHeaderV1 *header = new AliFastJetHeaderV1; | |
37 | header->SetBGMode(0); | |
38 | header->SetRadius(0.4); | |
39 | //header->SetGhostEtaMax(2); | |
40 | //header->SetGhostArea(0.05); | |
41 | header->SetAlgorithm(2); // antikt_algorithm = 2, kt = 0 (see fastjet/fastjet/JetDefinition.hh | |
42 | ||
43 | AliFastJetFinder *FastJet = new AliFastJetFinder; | |
44 | FastJet->SetJetHeader(header); | |
45 | ||
46 | AliAODEvent *aod = new AliAODEvent(); | |
47 | aod->CreateStdContent(); | |
48 | FastJet->ConnectAOD(aod); | |
49 | ||
50 | AliJetCalTrkEvent JetFinderEvent(0,1); | |
51 | TClonesArray *plist = new TClonesArray("TParticle"); | |
027fca2e | 52 | TClonesArray aliplist("AliMCParticle",1000); |
c45ff98f | 53 | |
54 | for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) { | |
55 | pythia->GenerateEvent(); | |
56 | ||
57 | pythia->GetParticles(plist); | |
58 | ||
027fca2e | 59 | aliplist.Clear(); |
60 | JetFinderEvent.Clear(); | |
c45ff98f | 61 | |
62 | Int_t n_part = plist->GetEntries(); | |
63 | for (Int_t i_part = 0; i_part < n_part; i_part++) { | |
64 | part=(TParticle*) plist->At(i_part); | |
65 | ||
66 | if (part->GetStatusCode() >= 10) // Not a final state particle | |
67 | continue; | |
68 | ||
027fca2e | 69 | if (i_part >= aliplist.GetSize()) { |
70 | aliplist.Expand(1.5*aliplist.GetSize()); | |
71 | } | |
c45ff98f | 72 | new (aliplist[i_part]) AliMCParticle(part); |
73 | ||
74 | JetFinderEvent.AddCalTrkTrackKine((AliMCParticle*)aliplist[i_part],1,1); | |
75 | } | |
76 | ||
77 | aod->ClearStd(); | |
78 | FastJet->Reset(); | |
79 | FastJet->SetCalTrkEvent(JetFinderEvent); | |
80 | FastJet->ProcessEvent(); | |
81 | if (aod->GetNJets() > 0) { | |
027fca2e | 82 | cout << "event " << iEvent << " " << aod->GetNJets() << " jets found" << endl; |
c45ff98f | 83 | for (Int_t iJet = 0; iJet < aod->GetNJets(); iJet++) { |
84 | jet = aod->GetJet(iJet); | |
85 | cout << "\t jet " << iJet << " pt " << jet->Pt() << " eta " << jet->Eta() << " phi " << jet->Phi() << endl; | |
86 | } | |
87 | } | |
88 | } | |
89 | } |