Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / JETAN / read_jets.C
1
2 void 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     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);
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));
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);
117                   e4H->Fill(emax);
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;
141
142         Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
143         
144         if (r  < rmin) {
145                 rmin = r;
146                 igen = j;
147         }
148     }
149     if (egen > 125. && egen < 150.) 
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");
174   dr2H->Draw();
175   dr3H->Draw("same");
176
177   TCanvas* c4 = new TCanvas("c4");
178   eH->Draw();
179
180   TCanvas* c5 = new TCanvas("c5");
181   etaH->Draw();
182
183   TCanvas* c5a = new TCanvas("c5a");
184   eta1H->Draw();
185
186   TCanvas* c5b = new TCanvas("c5b");
187   eta2H->Draw();
188
189   TCanvas* c6 = new TCanvas("c6");
190   e4H->Draw();
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();
198 }