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 }