Temporary removing AliEMCALv3 and AliEMCALHitv1 from the compilation until they are...
[u/mrichter/AliRoot.git] / JETAN / read_jets.C
CommitLineData
99e5fe42 1
2void read_jets(const char* fn = "jets.root")
3
4{
5 //
6 // Some histos
7 //
8 TH1F* eH = new TH1F("eH" , "Jet Energy", 40., 0., 200.);
9 TH1F* e1H = new TH1F("e1H" , "Jet Energy", 40., 0., 200.);
10 TH1F* e2H = new TH1F("e2H" , "Jet Energy", 40., 0., 200.);
11 TH1F* e3H = new TH1F("e3H" , "Jet Energy", 40., 0., 200.);
12 TH1F* e4H = new TH1F("e4H" , "Jet Energy", 40., 0., 200.);
13 TH1F* dr1H = new TH1F("dr1H", "delta R", 160., 0., 2.);
14 TH1F* dr2H = new TH1F("dr2H", "delta R", 160., 0., 2.);
15 TH1F* dr3H = new TH1F("dr4H", "delta R", 160., 0., 2.);
16 TH1F* etaH = new TH1F("etaH", "eta", 160., -2., 2.);
17
18
19 // load jet library
20 gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libJETAN");
21
22 // open file
23 TFile *jFile = new TFile(fn);
24
25 // get jet header and display parameters
26 AliUA1JetHeader* jHeader =
27 (AliUA1JetHeader*) (jFile->Get("AliUA1JetHeader"));
28 jHeader->PrintParameters();
29
30 // get reader header and events to be looped over
31 AliJetESDReaderHeader *jReaderH =
32 (AliJetESDReaderHeader*)(jFile->Get("AliJetESDReaderHeader"));
33 Int_t first = jReaderH->GetFirstEvent();
34 Int_t last = jReaderH->GetLastEvent();
35 cout << "First event = " << first << " Last event = " << last << endl;
36
37
38 // loop over events
39 AliJet *jets, *gjets;
40 AliLeading *leading;
41
42 for (Int_t i=first; i< last; i++) {
43 cout << " Analyzing event " << i << endl;
44 // get next tree with AliJet
45 char nameT[100];
46 sprintf(nameT, "TreeJ%d",i);
47 TTree *jetT =(TTree *)(jFile->Get(nameT));
48 jetT->SetBranchAddress("FoundJet", &jets);
49 jetT->SetBranchAddress("GenJet", &gjets);
50 jetT->SetBranchAddress("LeadingPart", &leading);
51 jetT->GetEntry(0);
52
53//
54// Find the jet with the highest E_T
55//
56 Int_t njets = jets->GetNJets();
57
58 Float_t emax = 0.;
59 Int_t imax = -1;
60 for (Int_t j = 0; j < njets; j++) {
61 if (jets->GetPt(j) > emax && TMath::Abs(jets->GetEta(j)) < 0.5) {
62 emax = jets->GetPt(j);
63 imax = j;
64 }
65 }
66
67 if (imax == -1) {
68 e2H->Fill(gjets->GetPt(0));
69 } else {
70 eH->Fill(jets->GetPt(imax));
71 dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0));
72//
73// Find the generated jet closest to the reconstructed
74//
75
76 Float_t rmin;
77 Int_t igen;
78 Float_t etaj = jets->GetEta(imax);
79 Float_t phij = jets->GetPhi(imax);
80
81 Int_t ngen = gjets->GetNJets();
82 if (ngen != 0) {
83 rmin = 1.e6;
84 igen = -1;
85 for (Int_t j = 0; j < ngen; j++) {
86 Float_t etag = gjets->GetEta(j);
87 Float_t phig = gjets->GetPhi(j);
88 Float_t deta = etag - etaj;
89 Float_t dphi = TMath::Abs(phig - phij);
90 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
91 Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
92 if (r < rmin) {
93 rmin = r;
94 igen = j;
95 }
96 }
97
98 Float_t egen = gjets->GetPt(igen);
99 e1H->Fill(gjets->GetPt(igen));
100
101 if (egen > 105. && egen < 125.) {
102 e4H->Fill(emax);
103 if (rmin > 0.3) etaH->Fill(etaj);
104 dr2H->Fill(rmin);
105 }
106 }
107 }
108
109
110//
111// Leading particle
112//
113 Float_t etal = leading->GetLeading()->Eta();
114 Float_t phil = leading->GetLeading()->Phi();
115 Float_t el = leading->GetLeading()->E();
116 e3H->Fill(el);
117
118 Float_t rmin = 1.e6;
119 Int_t igen = -1;
120 Int_t ngen = gjets->GetNJets();
121 for (Int_t j = 0; j < ngen; j++) {
122 Float_t etag = gjets->GetEta(j);
123 Float_t phig = gjets->GetPhi(j);
124 Float_t deta = etag-etal;
125 Float_t dphi = TMath::Abs(phig - phil);
126 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
127 Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
128
129 if (r < rmin) {
130 rmin = r;
131 igen = j;
132 }
133 }
134
135 dr3H->Fill(rmin);
136
137// cout << " Generated Jets:" << endl;
138// gjets->PrintJets();
139// cout << " Leading particle: " << endl;
140// leading->PrintLeading();
141 }
142
143 // Get Control Plots
144 TCanvas* c1 = new TCanvas("c1");
145 eH->Draw();
146 e1H->SetLineColor(2);
147 e2H->SetLineColor(4);
148 e3H->SetLineColor(5);
149 e1H->Draw("same");
150 e2H->Draw("same");
151 e3H->Draw("same");
152
153 TCanvas* c2 = new TCanvas("c2");
154// dr1H->Draw();
155 dr2H->SetLineColor(2);
156 dr2H->Draw("");
157
158 TCanvas* c3 = new TCanvas("c3");
159 dr3H->Draw();
160 dr2H->Draw("same");
161
162 TCanvas* c4 = new TCanvas("c4");
163 eH->Draw();
164
165 TCanvas* c5 = new TCanvas("c5");
166 etaH->Draw();
167
168 TCanvas* c6 = new TCanvas("c6");
169 e4H->Draw();
170}