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 }