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 }