EMCAL/AnaJets.C

   1 void AnaJets(Int_t evNumber1=0, Int_t evNumber2=0)
   2 {
   3 //*-- Author: Andreas Morsch (CERN)
   4
   5     if (gClassTable->GetID("AliRun") < 0) {
   6         gROOT->LoadMacro("loadlibs.C");
   7         loadlibs();
   8     }
   9 // Connect the Root Galice file containing Geometry, Kine and Hits
  10
  11     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  12
  13     if (!file) {
  14         printf("\n Creating galice.root \n");
  15         file = new TFile("galice.root");
  16     } else {
  17         printf("\n galice.root found in file list");
  18     }
  19 // Get AliRun object from file or create it if not on file
  20     if (!gAlice) {
  21         gAlice = (AliRun*)(file->Get("gAlice"));
  22         if (gAlice) printf("AliRun object found on file\n");
  23         if (!gAlice) {
  24             printf("\n create new gAlice object");
  25             gAlice = new AliRun("gAlice","Alice test program");
  26         }
  27     }
  28 // Book histos
  29     TH1F *eH   = new TH1F("eH","Energy",    150,  0.0, 150.);
  30     TH1F *etaH = new TH1F("eEta","Eta",     180, -0.9,  0.9);
  31     TH1F *phiH = new TH1F("ePhi","Phi",      62, -3.1,  3.1);
  32     TH1F *tH   = new TH1F("tH","n tracks",   30,  0.5, 29.5);
  33     TH1F *ptH  = new TH1F("ptH","Track pT", 100., 0., 100.);
  34     TH1F *drH  = new TH1F("drH","Track dR", 120., 0.,   6.);
  35
  36     Float_t phiT[50], etaT[50], ptT[50];
  37
  38
  39     TClonesArray* jets = new TClonesArray("AliEMCALJet",10000);
  40
  41     for (int nev=0; nev<= evNumber2; nev++) {
  42         printf("\n Event .............%d", nev);
  43         Int_t nparticles = gAlice->GetEvent(nev);
  44         Int_t nbytes     = 0;
  45         AliEMCAL *pEMCAL  = (AliEMCAL*) gAlice->GetModule("EMCAL");
  46         if (pEMCAL) {
  47             TTree *TR = gAlice->TreeR();
  48             Int_t nent=TR->GetEntries();
  49             TR->SetBranchAddress("EMCALJets", &jets);
  50             nbytes += TR->GetEntry(0);
  51             Int_t nJet = jets->GetEntries();
  52             printf("\n Number of Jets %d", nJet);
  53             AliEMCALJet  *mJet;
  54             for (Int_t ij=0; ij < nJet; ij++) {
  55                 mJet = (AliEMCALJet*)jets->UncheckedAt(ij);
  56                 Float_t eta = mJet->Eta();
  57
  58                 printf("\n Jet:%d E %f phi %f eta %f tracks %d\n", ij,
  59                        mJet->Energy(), mJet->Phi(), eta,
  60                        mJet->NTracks());
  61                 etaH->Fill(mJet->Eta());
  62                 phiH->Fill(mJet->Phi());
  63                 if (TMath::Abs(eta) < 0.4)
  64                     eH->Fill(mJet->Energy());
  65                 tH  ->Fill((Float_t)mJet->NTracks());
  66
  67                 mJet->TrackList(ptT, etaT, phiT);
  68                 for (Int_t it = 0; it < mJet->NTracks(); it++)
  69                 {
  70                     printf("\n Track: %5d pT %8.3f eta %8.3f phi %8.3f",
  71                            it, ptT[it], etaT[it], phiT[it]);
  72                     ptH->Fill(ptT[it]);
  73                     Float_t dPhi = phiT[it]-mJet->Phi();
  74                     Float_t dEta = etaT[it]-mJet->Eta();
  75                     Float_t dr = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
  76                     drH->Fill(dr);
  77                 }
  78            } // jet
  79        } // ?EMCAL
  80    } // event
  81     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
  82     c1->Divide(2,2);
  83     c1->cd(1); eH->Draw();
  84     c1->cd(2); etaH->Draw();
  85     c1->cd(3); phiH->Draw();
  86     c1->cd(4); tH->Draw();
  87
  88     TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
  89     c2->Divide(2,2);
  90     c2->cd(1); ptH->Draw();
  91     c2->cd(2); drH->Draw();
  92 }
  93