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