First commit of new jet reconstruction and analysis code to be used for the
[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     
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 }