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