correcting baryon mass subtraction for visible energy in MC
[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.);
83a444b1 16 TH1F* etaH = new TH1F("etaH", "eta", 160., -2., 2.);
17 TH1F* eta1H = new TH1F("eta1H", "eta", 160., -2., 2.);
18 TH1F* eta2H = new TH1F("eta2H", "eta", 160., -2., 2.);
19
20 TH1F* phiH = new TH1F("phiH", "phi", 160., -3., 3.);
21 TH1F* phi1H = new TH1F("phi1H", "phi", 160., 0., 6.28);
22 TH1F* phi2H = new TH1F("phi2H", "phi", 160., 0., 6.28);
99e5fe42 23
24
25 // load jet library
26 gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libJETAN");
27
28 // open file
29 TFile *jFile = new TFile(fn);
30
31 // get jet header and display parameters
32 AliUA1JetHeader* jHeader =
33 (AliUA1JetHeader*) (jFile->Get("AliUA1JetHeader"));
34 jHeader->PrintParameters();
35
36 // get reader header and events to be looped over
37 AliJetESDReaderHeader *jReaderH =
38 (AliJetESDReaderHeader*)(jFile->Get("AliJetESDReaderHeader"));
39 Int_t first = jReaderH->GetFirstEvent();
40 Int_t last = jReaderH->GetLastEvent();
41 cout << "First event = " << first << " Last event = " << last << endl;
42
43
44 // loop over events
45 AliJet *jets, *gjets;
46 AliLeading *leading;
47
48 for (Int_t i=first; i< last; i++) {
49 cout << " Analyzing event " << i << endl;
50 // get next tree with AliJet
51 char nameT[100];
52 sprintf(nameT, "TreeJ%d",i);
53 TTree *jetT =(TTree *)(jFile->Get(nameT));
54 jetT->SetBranchAddress("FoundJet", &jets);
55 jetT->SetBranchAddress("GenJet", &gjets);
56 jetT->SetBranchAddress("LeadingPart", &leading);
57 jetT->GetEntry(0);
58
59//
60// Find the jet with the highest E_T
61//
62 Int_t njets = jets->GetNJets();
63
64 Float_t emax = 0.;
65 Int_t imax = -1;
66 for (Int_t j = 0; j < njets; j++) {
67 if (jets->GetPt(j) > emax && TMath::Abs(jets->GetEta(j)) < 0.5) {
68 emax = jets->GetPt(j);
69 imax = j;
70 }
71 }
72
73 if (imax == -1) {
74 e2H->Fill(gjets->GetPt(0));
75 } else {
76 eH->Fill(jets->GetPt(imax));
77 dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0));
78//
79// Find the generated jet closest to the reconstructed
80//
81
82 Float_t rmin;
83 Int_t igen;
84 Float_t etaj = jets->GetEta(imax);
85 Float_t phij = jets->GetPhi(imax);
86
87 Int_t ngen = gjets->GetNJets();
88 if (ngen != 0) {
89 rmin = 1.e6;
90 igen = -1;
91 for (Int_t j = 0; j < ngen; j++) {
92 Float_t etag = gjets->GetEta(j);
93 Float_t phig = gjets->GetPhi(j);
94 Float_t deta = etag - etaj;
95 Float_t dphi = TMath::Abs(phig - phij);
96 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
97 Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
98 if (r < rmin) {
99 rmin = r;
100 igen = j;
101 }
102 }
103
104 Float_t egen = gjets->GetPt(igen);
105 e1H->Fill(gjets->GetPt(igen));
83a444b1 106 Float_t etag = gjets->GetEta(igen);
107 Float_t phig = gjets->GetPhi(igen);
108 Float_t dphi = phig - phij;
109
110 if (egen > 125. && egen < 150.) {
111 phiH->Fill((dphi));
112 etaH->Fill(etag - etaj);
113 phi1H->Fill(phig);
114 phi2H->Fill(phij);
115 eta1H->Fill(etag);
116 eta2H->Fill(etaj);
99e5fe42 117 e4H->Fill(emax);
99e5fe42 118 dr2H->Fill(rmin);
119 }
120 }
121 }
122
123
124//
125// Leading particle
126//
127 Float_t etal = leading->GetLeading()->Eta();
128 Float_t phil = leading->GetLeading()->Phi();
129 Float_t el = leading->GetLeading()->E();
130 e3H->Fill(el);
131
132 Float_t rmin = 1.e6;
133 Int_t igen = -1;
134 Int_t ngen = gjets->GetNJets();
135 for (Int_t j = 0; j < ngen; j++) {
136 Float_t etag = gjets->GetEta(j);
137 Float_t phig = gjets->GetPhi(j);
138 Float_t deta = etag-etal;
139 Float_t dphi = TMath::Abs(phig - phil);
140 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
83a444b1 141
99e5fe42 142 Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
143
144 if (r < rmin) {
145 rmin = r;
146 igen = j;
147 }
148 }
83a444b1 149 if (egen > 125. && egen < 150.)
99e5fe42 150 dr3H->Fill(rmin);
151
152// cout << " Generated Jets:" << endl;
153// gjets->PrintJets();
154// cout << " Leading particle: " << endl;
155// leading->PrintLeading();
156 }
157
158 // Get Control Plots
159 TCanvas* c1 = new TCanvas("c1");
160 eH->Draw();
161 e1H->SetLineColor(2);
162 e2H->SetLineColor(4);
163 e3H->SetLineColor(5);
164 e1H->Draw("same");
165 e2H->Draw("same");
166 e3H->Draw("same");
167
168 TCanvas* c2 = new TCanvas("c2");
169// dr1H->Draw();
170 dr2H->SetLineColor(2);
171 dr2H->Draw("");
172
173 TCanvas* c3 = new TCanvas("c3");
83a444b1 174 dr2H->Draw();
175 dr3H->Draw("same");
99e5fe42 176
177 TCanvas* c4 = new TCanvas("c4");
178 eH->Draw();
179
180 TCanvas* c5 = new TCanvas("c5");
181 etaH->Draw();
182
83a444b1 183 TCanvas* c5a = new TCanvas("c5a");
184 eta1H->Draw();
185
186 TCanvas* c5b = new TCanvas("c5b");
187 eta2H->Draw();
188
99e5fe42 189 TCanvas* c6 = new TCanvas("c6");
190 e4H->Draw();
83a444b1 191 TCanvas* c7 = new TCanvas("c7");
192 phiH->Draw();
193
194 TCanvas* c7a = new TCanvas("c7a");
195 phi1H->Draw();
196 TCanvas* c7b = new TCanvas("c7b");
197 phi2H->Draw();
99e5fe42 198}