EMCAL/AliEMCALJetMicroDst.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 /* $Id$ */
  18
  19 //*-- Authors: Aleksei Pavlinov (WSU)
  20 //*
  21
  22 #include "AliEMCALJetMicroDst.h"
  23 #include "AliRun.h"
  24 #include "AliHeader.h"
  25 #include "AliGenEventHeader.h"
  26 #include "AliGenHijingEventHeader.h"
  27 #include "AliEMCALJetFinder.h"
  28 #include <TVector3.h>
  29 #include <TROOT.h>
  30 #include <TBrowser.h>
  31 #include <TString.h>
  32 #include <TParticle.h>
  33
  34 ClassImp(AliEMCALJetMicroDst)
  35
  36 TString nameTree("jetMDST"); // 7-feb-2002
  37
  38 TH1F*  hPtPart, *hNJet, *hPtJet;
  39 TH2F*  hEtaPhiPart, *hEtaPhiJet;
  40
  41 AliEMCALJetMicroDst::AliEMCALJetMicroDst(char *name, char *tit) : TNamed(name,tit)
  42 {
  43   fFile  = 0;
  44   fTree  = 0;
  45   fDebug = 0;
  46   fName  = "dstTest.root"; // default name
  47
  48 //  Don't add histos to the current directory
  49 //  TH1::AddDirectory(0);
  50   gROOT->cd();
  51   hPtPart     = new TH1F("hPtPart","P_{T} for partons", 300, 0., 300.);
  52   hEtaPhiPart = new TH2F("hEtaPhiPart","#eta #phi distr.for partons after HSc",
  53                        28, -0.7, 0.7, 21, 0.0, TMath::Pi()*2./3.);
  54
  55   hNJet  = new TH1F("hNJet","number of jets", 11, -0.5, 10.5);
  56   hPtJet = new TH1F("hPtJet","P_{T} for jets", 300, 0., 300.);
  57   hEtaPhiJet = new TH2F("hEtaPhiJet","#eta #phi distr.for jets (W)",
  58                        28, -0.7, 0.7, 21, 0.0, TMath::Pi()*2./3.);
  59   fListHist = MoveHistsToList("Hist For AliEMCALJetMicroDst", kFALSE);
  60 }
  61
  62 AliEMCALJetMicroDst::~AliEMCALJetMicroDst()
  63 {
  64   if(fFile) fFile->Close();
  65 }
  66
  67 Bool_t AliEMCALJetMicroDst::Create(TFile *file)
  68 {
  69   if(!file) {
  70     printf("<E> AliEMCALJetMicroDst::create -> define TFile for output\n");
  71     return kFALSE;
  72   }
  73   fFile = file;
  74   fFile->cd();
  75   fTree = new TTree(nameTree.Data(),"Temporary micro DST for jets analysis");
  76   // partons
  77   fTree->Branch("npart", &fNpart, "npart/I");
  78   fTree->Branch("xpt",  fXpt,  "xpt[npart]/F");
  79   fTree->Branch("xeta", fXeta, "xeta[npart]/F");
  80   fTree->Branch("xphi", fXphi, "xphi[npart]/F");
  81   // jets
  82   fTree->Branch("njet", &fNjet, "njet/I");
  83   fTree->Branch("jet",  fJet,  "jet[njet]/F");
  84   fTree->Branch("jetaw", fJetaW, "jetaw[njet]/F");
  85   fTree->Branch("jphiw", fJphiW, "jphiw[njet]/F");
  86   fTree->Branch("jetal", fJetaL, "jetal[njet]/F");
  87   fTree->Branch("jphil", fJphiL, "jphil[njet]/F");
  88   return kTRUE;
  89 }
  90
  91 Bool_t AliEMCALJetMicroDst::Create(const char *fname)
  92 {
  93   TFile *file = new TFile(fname, "RECREATE");
  94   if(file) {
  95     fName = fname;
  96     return Create(file);
  97   } else     return kFALSE;
  98 }
  99
 100 Bool_t AliEMCALJetMicroDst::Open(const char *fname)
 101 {
 102   if(fFile && fFile->IsOpen()) fFile->Close();
 103   if(strlen(fname)) fName = fname;
 104   TFile *file = new TFile(fName.Data(), "READ");
 105   if(file) {
 106     printf("<I> open file %s \n",fName.Data());
 107     return Initialize(file);
 108   } else {
 109     printf("<E> can not open file %s \n",fName.Data());
 110     return kFALSE;
 111   }
 112 }
 113
 114 Bool_t AliEMCALJetMicroDst::Initialize(TFile *file)
 115 {
 116   if(file) fFile = file;
 117   fFile->cd();
 118   fTree = (TTree*)fFile->Get(nameTree.Data());
 119   if(!fTree) return kFALSE;
 120   // partons
 121   fTree->SetBranchAddress("npart",&fNpart);
 122   fTree->SetBranchAddress("xpt",   fXpt);
 123   fTree->SetBranchAddress("xeta",  fXeta);
 124   fTree->SetBranchAddress("xphi",  fXphi);
 125   // jets
 126   fTree->SetBranchAddress("njet", &fNjet);
 127   fTree->SetBranchAddress("jet",   fJet);
 128   fTree->SetBranchAddress("jetaw", fJetaW);
 129   fTree->SetBranchAddress("jphiw", fJphiW);
 130   fTree->SetBranchAddress("jetal", fJetaL);
 131   fTree->SetBranchAddress("jphil", fJphiL);
 132
 133   return kTRUE;
 134 }
 135
 136 void AliEMCALJetMicroDst::Print(Option_t* option) const
 137 {
 138   if(option);
 139   if(fFile) {
 140     fFile->Print();
 141     if(fTree) fTree->Print();
 142     else printf("<I> TRee is zero\n");
 143   } else printf("<I> File with TRee is zero\n");
 144   printf("Default name of file %s\n", fName.Data());
 145
 146   printf("\n #partons %2i ", fNpart);
 147   for(Int_t i=0; i<fNpart; i++){
 148     printf("\n     %1i Pt %7.1f eta  %f7.4 phi  %f7.4 ",
 149     i, fXpt[i],fXeta[i],fXphi[i]);
 150   }
 151
 152   printf("\n #jets    %2i ", fNjet);
 153   for(Int_t i=0; i<fNjet; i++){
 154     printf("\n     %1i Et %7.1f etaW %f7.4 phiW %f7.4 ",
 155     i,fJet[i],fJetaW[i],fJphiW[i]);
 156   }
 157 }
 158
 159 void AliEMCALJetMicroDst::Fill(AliRun *run, AliEMCALJetFinder* jetFinder)
 160 {
 161   if(!run) run = gAlice;
 162   AliGenEventHeader* evHeader = run->GetHeader()->GenEventHeader();
 163   TString tmp(evHeader->ClassName());
 164   if(tmp.Contains("Hijing")) {
 165      AliGenHijingEventHeader *hijEvHeader = (AliGenHijingEventHeader*)evHeader;
 166      FillPartons(hijEvHeader);
 167   } else if(tmp.Contains("Pythia")) {
 168      FillPartons();
 169   } else {
 170     printf("\n <E> Wrong type of generator -> %s ",tmp.Data());
 171     printf("\n     Inof about partons will be absent ");
 172   }
 173
 174   FillJets(jetFinder);
 175   fTree->Fill();
 176 }
 177
 178 void AliEMCALJetMicroDst::FillPartons(AliGenHijingEventHeader *header)
 179 {
 180   TLorentzVector parton[4];
 181   header->GetJets(parton[0], parton[1], parton[2], parton[3]);
 182
 183   fNpart = 4; //
 184   for(Int_t i=0; i<4; i++){
 185      fXpt[i]  = parton[i].Pt();
 186      fXeta[i] = parton[i].Eta();
 187      fXphi[i] = parton[i].Phi();
 188      if(i>1) {
 189        hEtaPhiPart->Fill(parton[i].Eta(), parton[i].Phi());
 190        hPtPart->Fill(parton[i].Pt());
 191      }
 192   }
 193 }
 194
 195 void AliEMCALJetMicroDst::FillPartons()
 196 {// for case of Pythia -> get info from full event record
 197
 198   fNpart = 2;
 199   TParticle *MPart;
 200   Int_t ind;
 201   for(Int_t i=6; i<8; i++){
 202      MPart = gAlice->Particle(i);
 203      ind   = i-6;
 204      fXpt[ind]  = MPart->Pt();
 205      fXeta[ind] = MPart->Eta();
 206      fXphi[ind] = MPart->Phi();
 207
 208      hEtaPhiPart->Fill(fXeta[ind], fXphi[ind]);
 209      hPtPart->Fill(fXpt[ind]);
 210   }
 211 }
 212
 213 void AliEMCALJetMicroDst::FillJets(AliEMCALJetFinder* jetFinder)
 214 {
 215   fNjet = 0;
 216   if(fDebug>1) printf("\n<I> AliEMCALJetMicroDst::FillJets");
 217   if(!jetFinder) {
 218     if(fDebug>1) printf("\n : jetFinder is zero");
 219     return;
 220   }
 221   fNjet = jetFinder->Njets();
 222   if(fNjet>10) {
 223     if(fDebug>1) printf("\n <W> wrong value of jetFinder->Njets() %i ", fNjet);
 224     fNjet = 10;
 225   }
 226   hNJet->Fill(fNjet);
 227   if(fDebug>1) printf("\n <I> fNjet %i", fNjet);
 228   if(fNjet){
 229     for(Int_t i=0; i<fNjet; i++){
 230       fJet[i]   = jetFinder->JetEnergy(i);
 231       fJetaW[i] = jetFinder->JetEtaW(i);
 232       fJphiW[i] = jetFinder->JetPhiW(i);
 233       fJetaL[i] = jetFinder->JetEtaL(i);
 234       fJphiL[i] = jetFinder->JetPhiL(i);
 235       hPtJet->Fill(fJet[i]);
 236       hEtaPhiJet->Fill(fJetaW[i],fJphiW[i]);
 237     }
 238   }
 239 }
 240
 241 Int_t AliEMCALJetMicroDst::GetEntry(Int_t entry)
 242 { // Read contents of entry.
 243    if (!fTree) {
 244       printf("\n<E> AliEMCALJetMicroDst::GetEntry() -> define TTree \n");
 245       return -1;
 246    }
 247    return fTree->GetEntry(entry);
 248 }
 249
 250 Bool_t AliEMCALJetMicroDst::GetParton(Int_t i, Float_t& pt, Float_t& eta, Float_t& phi)
 251 {
 252   if(i>=0 && i<fNpart) {
 253     pt  = fXpt[i];
 254     eta = fXeta[i];
 255     phi = fXphi[i];
 256     return kTRUE;
 257   } else return kFALSE;
 258 }
 259
 260 Bool_t AliEMCALJetMicroDst::GetParton(Int_t i, TVector3& vec)
 261 {
 262   static Float_t pt, eta, phi;
 263
 264   if(!GetParton(i, pt, eta, phi)) return kFALSE;
 265
 266   FillVector(pt, eta, phi, vec);
 267   return kTRUE;
 268 }
 269
 270 Bool_t AliEMCALJetMicroDst::GetJet(Int_t i,Int_t mode,Float_t& pt, Float_t& eta, Float_t& phi)
 271 {// mode=1(W) mode=any(L)
 272   if(i>=0 && i<fNjet) {
 273     pt  = fJet[i];
 274     if(mode==1) {
 275       eta = fJetaW[i];
 276       phi = fJphiW[i];
 277     } else {
 278       eta = fJetaL[i];
 279       phi = fJphiL[i];
 280     }
 281     return kTRUE;
 282   } else return kFALSE;
 283 }
 284
 285 Bool_t AliEMCALJetMicroDst::GetJet(Int_t i, Int_t mode, TVector3& vec)
 286 {
 287   static Float_t pt, eta, phi;
 288
 289   if(!GetJet(i, mode, pt, eta, phi)) return kFALSE;
 290   FillVector(pt, eta, phi, vec);
 291   return kTRUE;
 292 }
 293
 294 void AliEMCALJetMicroDst::Test()
 295 {
 296   if(!fFile || !fTree ) {
 297     printf("\n<I> AliEMCALJetMicroDst::Test() -> define file with proper TTree !");
 298     return;
 299   }
 300   Int_t nbytes=0, nb=0, nentries=Int_t(fTree->GetEntries());
 301   for(Int_t i=0; i<nentries; i++){
 302     nb = fTree->GetEntry(i);
 303     nbytes += nb;
 304     for(Int_t j=0; j<fNpart; j++){
 305       hEtaPhiPart->Fill(fXeta[j], fXphi[j]);
 306       hPtPart->Fill(fXpt[j]);
 307     }
 308
 309     hNJet->Fill(fNjet);
 310     if(fNjet){
 311       for(Int_t j=0; j<fNjet; j++) {
 312         hPtJet->Fill(fJet[j]);
 313         hEtaPhiJet->Fill(fJetaW[j],fJphiW[j]);
 314       }
 315     }
 316   }
 317   printf("\n<I> AliEMCALJetMicroDst::Test() -> Entries %5i Bytes %10i\n", nentries, nb);
 318 }
 319
 320 void AliEMCALJetMicroDst::FillVector(Float_t pt, Float_t eta, Float_t phi, TVector3& vec)
 321 { // service function
 322   static Float_t px, py, pz;
 323
 324   px = pt*TMath::Cos(phi);
 325   py = pt*TMath::Sin(phi);
 326   pz = pt*TMath::SinH(eta);  // sinh(eta) = cot(theta)
 327
 328   vec.SetXYZ(px, py, pz);
 329 }
 330
 331 void AliEMCALJetMicroDst::Close()
 332 {
 333   fFile->Write();
 334   fTree->Print();
 335   fFile->Close();
 336   fFile = 0;
 337   fTree = 0;
 338 }
 339
 340 void AliEMCALJetMicroDst::Browse(TBrowser* b)
 341 {
 342    if(fTree)      b->Add((TObject*)fTree);
 343    if(fListHist)  b->Add((TObject*)fListHist);
 344    //   TObject::Browse(b);
 345 }
 346
 347 Bool_t  AliEMCALJetMicroDst::IsFolder() const
 348 {
 349   if(fTree || fListHist) return kTRUE;
 350   else                   return kFALSE;
 351 }
 352
 353 TList* AliEMCALJetMicroDst::MoveHistsToList(char* name, Bool_t putToBrowser)
 354 {
 355   gROOT->cd();
 356   TIter nextHist(gDirectory->GetList());
 357   TList *list = new TList;
 358   list->SetName(name);
 359   TObject *objHist;
 360   while((objHist=nextHist())){
 361     if (!objHist->InheritsFrom("TH1")) continue;
 362     ((TH1*)objHist)->SetDirectory(0); // Remove from gROOT
 363     list->Add(objHist);
 364   }
 365   if(putToBrowser) gROOT->GetListOfBrowsables()->Add((TObject*)list);
 366   return list;
 367 }