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