Updating example reco macro for real data
[u/mrichter/AliRoot.git] / JETAN / read_jets_kine.C
1
2 void 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 }