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
  23 #include "AliEMCALJetMicroDst.h"
  24 #include "AliRun.h"
  25 #include "AliHeader.h"
  26 #include "AliGenEventHeader.h"
  27 #include "AliGenHijingEventHeader.h"
  28 #include "AliEMCALJetFinder.h"
  29 #include <TVector3.h>
  30 #include <TROOT.h>
  31 #include <TBrowser.h>
  32 #include <TString.h>
  33 #include <TParticle.h>
  34 #include "AliMC.h"
  35 #include <TH2.h>
  36 #include <TH1F.h>
  37
  38
  39 ClassImp(AliEMCALJetMicroDst)
  40
  41 TString gAliNameTree("jetMDST"); // 7-feb-2002
  42
  43 //TH1F*  fhPtPart, *fhNJet, *fhPtJet;
  44 //TH2F*  fhEtaPhiPart, *fhEtaPhiJet;
  45 //TH1F*  fhNcell, *fhCellId, *fhCellEt, *fhSumEt;
  46 //TH1F*  fhNgrid, *fhGridId, *fhGridEt, *fhSumEtGrForJF;
  47
  48 extern "C" void sgpdge_(Int_t &i, Int_t &pdggea);
  49
  50 AliEMCALJetMicroDst::AliEMCALJetMicroDst(const char *name, const char *tit) : TNamed(name,tit)
  51 {
  52         //constructor
  53   fFile  = 0;
  54   fTree  = 0;
  55   fDebug = 0;
  56
  57 //  Don't add histos to the current directory
  58 //  TH1::AddDirectory(0);
  59   gROOT->cd();
  60   fhPtPart     = new TH1F("fhPtPart","P_{T} for partons", 300, 0., 300.);
  61   // 16-jan-2002 - new limit fo phi
  62   fhEtaPhiPart = new TH2F("fhEtaPhiPart","#eta #phi distr.for partons after HSc",
  63                        28, -0.7, 0.7, 21, 0.0, (2.0/3.0)*TMath::Pi());
  64
  65   fhNJet  = new TH1F("fhNJet","number of jets", 11, -0.5, 10.5);
  66   fhPtJet = new TH1F("fhPtJet","P_{T} for jets", 500, 0., 500.);
  67   fhEtaPhiJet = new TH2F("fhEtaPhiJet","#eta #phi distr.for jets (W)",
  68                        28, -0.7, 0.7, 21, 0.0, (2.0/3.0)*TMath::Pi());
  69
  70   fhNcell  = new TH1F("fhNcell","#cell with de>0.0 for EMCAL", 1400, 0.0, 14000.);
  71   fhCellId = new TH1F("fhCellId","cell ID with de>0.0 for EMCAL", 1400, 0.0, 14000.);
  72   fhCellEt = new TH1F("fhCellEt","cell Et for EMCAL", 1000, 0.0, 10.);
  73   fhSumEt  = new TH1F("fhSumEt","sum Et for EMCAL", 1000, 0.0, 1000.);
  74
  75   fhNgrid  = new TH1F("fhNgrid","#cell with de>0.0 in EMCAL grid for JF", 1400, 0.0, 14000.);
  76   fhGridId = new TH1F("fhGridId","cell ID with de>0.0 in EMCAL grid for JF", 1400, 0.0, 14000.);
  77   fhGridEt = new TH1F("fhGridEt","cell Et in EMCAL grid for JF", 1000, 0.0, 10.);
  78   fhSumEtGrForJF  = new TH1F("fhSumEtGrForJF","sum Et in EMCAL grid for JF", 1000, 0.0, 1000.);
  79
  80   fListHist = MoveHistsToList("Hist For AliEMCALJetMicroDst", kFALSE);
  81 }
  82
  83 AliEMCALJetMicroDst::~AliEMCALJetMicroDst()
  84 {
  85         //destructor
  86   if(fFile) fFile->Close();
  87 }
  88
  89 Bool_t AliEMCALJetMicroDst::Create(TFile *file)
  90 {
  91   // Creates the DST file
  92   if(!file) {
  93     Error("Create", "define TFile for output\n");
  94     return kFALSE;
  95   }
  96   fFile     = file;
  97   fFileName = fFile->GetName();
  98   fFile->cd();
  99   fTree = new TTree(gAliNameTree.Data(),"Temporary micro DST for jets analysis");
 100   // for jet calibration - 4-mar-2003
 101   fTree->Branch("fdecone", &fdecone, "fdecone/F");
 102   fTree->Branch("fptcone", &fptcone, "fptcone/F");
 103   // partons
 104   fTree->Branch("fnpart", &fnpart, "fnpart/I");
 105   fTree->Branch("fxpt",  fxpt,  "fxpt[fnpart]/F");
 106   fTree->Branch("fxeta", fxeta, "fxeta[fnpart]/F");
 107   fTree->Branch("fxphi", fxphi, "fxphi[fnpart]/F");
 108   // jets
 109   fTree->Branch("fnjet", &fnjet, "fnjet/I");
 110   fTree->Branch("fjet",  fjet,  "fjet[fnjet]/F");
 111   fTree->Branch("fjetaw", fjetaw, "fjetaw[fnjet]/F");
 112   fTree->Branch("fjphiw", fjphiw, "fjphiw[fnjet]/F");
 113   fTree->Branch("fjetal", fjetal, "fjetal[fnjet]/F");
 114   fTree->Branch("fjphil", fjphil, "fjphil[fnjet]/F");
 115
 116   // Et in EMCAL itself
 117   fTree->Branch("fncell", &fncell, "fncell/I");
 118   fTree->Branch("fidcell", fidcell, "fidcell[fncell]/I");
 119   fTree->Branch("fetcell", fetcell, "fetcell[fncell]/F");
 120
 121   // Et in EMCAL grid for JF
 122   fTree->Branch("fngrid", &fngrid,  "fngrid/I");
 123   fTree->Branch("fidgrid", fidgrid, "fidgrid[fngrid]/I");
 124   fTree->Branch("fetgrid", fetgrid, "fetgrid[fngrid]/F");
 125
 126   // charge particle which hit to EMCAL
 127   fTree->Branch("fnchp", &fnchp, "fnchp/I");
 128   fTree->Branch("fpid", fpid, "fpid[fnchp]/I");
 129   fTree->Branch("fppt", fppt, "fppt[fnchp]/F");
 130   fTree->Branch("fpeta", fpeta, "fpeta[fnchp]/F");
 131   fTree->Branch("fpphi", fpphi, "fpphi[fnchp]/F");
 132
 133   return kTRUE;
 134 }
 135
 136 Bool_t AliEMCALJetMicroDst::Create(const char *fname)
 137 {
 138         // Create member
 139   TFile *file = new TFile(fname, "RECREATE");
 140   if(file) {
 141     //    fNameFile = fname;
 142     return Create(file);
 143   } else   return kFALSE;
 144 }
 145
 146 Bool_t AliEMCALJetMicroDst::Open(const char *fname)
 147 {
 148         //Open member
 149   if(fFile && fFile->IsOpen()) fFile->Close();
 150   if(strlen(fname)) fName = fname;
 151   TFile *file = new TFile(fName.Data(), "READ");
 152   if(file) {
 153     Bool_t ini =  Initialize(file);
 154     Info("Open", "open file %s : initialize TTree %i",fName.Data(), Int_t(ini));
 155     return ini;
 156   } else {
 157     Error("Open", "can not open file %s",fName.Data());
 158     return kFALSE;
 159   }
 160 }
 161
 162 const Char_t* AliEMCALJetMicroDst::DefineName(Int_t mode)
 163 {
 164         //DefineName member
 165   static TString dir, name;
 166   //  dir = "jetDST/"; // 24-jan-2003
 167   dir = "/auto/alice/pavlinov/jet/microDST/"; // 24-jan-2003
 168   switch (mode) {
 169   case 1: // for characteristic of BG
 170     name = dir + "mDst1_1.root"; // Bg 2000 - first version
 171     SetTitle("Bg2000");
 172     break;
 173   case 2: // for characteristic of BG
 174     name = dir + "mDst2_1.root"; // Bg 4000 - first version
 175     SetTitle("Bg4000");
 176     break;
 177   case 3: // for characteristic of BG
 178     name = dir + "mDst3_1.root"; // Bg 8000 - first version
 179     SetTitle("Bg8000");
 180     break;
 181   case 4: // Central Hijing - 18-mar-2003
 182     name  = dir  + "march/";
 183     name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0HijingCentral.root";
 184     SetTitle("HijingCentral");
 185     break;
 186   case 5: // Para Hijing (Dn/Dy=8000) - 21-mar-2003
 187     name  = dir  + "march/";
 188     name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0ParaHijing8000.root";
 189     SetTitle("HIJINGparaDnDy8000");
 190     break;
 191   case 6: // Para Hijing (Dn/Dy=4000) - 21-mar-2003
 192     name  = dir  + "march/";
 193     name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0ParaHijing4000.root";
 194     SetTitle("HIJINGparaDnDy4000");
 195     break;
 196   case 7: // Para Hijing (Dn/Dy=2000) - 21-mar-2003
 197     name  = dir  + "march/";
 198     name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0ParaHijing2000.root";
 199     SetTitle("HIJINGparaDnDy2000");
 200     break;
 201   case 11: // pure PYTHIA with default value of parameters
 202     name = dir + "jF_R0.50MinCell0.0PtCut0.0EtSeed8.0MinEt40.0BGSubtr0SF11.6.root";
 203     break;
 204   case 12: // 0 + background
 205     name = dir + "jF_R0.50MinCell0.0PtCut0.0EtSeed8.0MinEt40.0BGSubtr0SF11.6kBackground2000.root";
 206     break;
 207     // calibration case
 208   case 101:
 209     name  = dir  + "march/";
 210     name += "Pythia100_1.root";
 211     SetTitle("Pythia100_1");
 212     break;
 213   case 102: // 2-apr-2003
 214     name  = dir  + "march/";
 215     name += "Pythia50_1.root";
 216     SetTitle("Pythia50_1");
 217     //    name = "microDst3th.root"; // 101 + (smearing and eff) - 14-mar-2003
 218     break;
 219   case 103:// 4-apr-2003
 220     name  = dir  + "march/";
 221     name += "Pythia200_1.root";
 222     SetTitle("Pythia200_1");
 223     //    name = "microDst4th.root"; // 102 + MinCell= 1.0
 224     break;
 225   default:
 226     Fatal("DefineName", "NO  D E F A U L T : mode %i\n", mode);
 227   }
 228   Info("DefineName", "mode %5i file : %s : Title %s\n", mode, name.Data(), GetTitle());
 229   return name.Data();
 230 }
 231
 232 Bool_t AliEMCALJetMicroDst::Initialize(TFile *file)
 233 {
 234         // Initialize method
 235   if(file) fFile = file;
 236   fFile->cd();
 237   fTree = (TTree*)fFile->Get(gAliNameTree.Data());
 238   if(!fTree) return kFALSE;
 239   // for jet calibration - 4-mar-2003
 240   fTree->SetBranchAddress("fdecone",&fdecone);
 241   fTree->SetBranchAddress("fptcone",&fptcone);
 242   // partons
 243   fTree->SetBranchAddress("fnpart",&fnpart);
 244   fTree->SetBranchAddress("fxpt",   fxpt);
 245   fTree->SetBranchAddress("fxeta",  fxeta);
 246   fTree->SetBranchAddress("fxphi",  fxphi);
 247   // jets
 248   fTree->SetBranchAddress("fnjet", &fnjet);
 249   fTree->SetBranchAddress("fjet",   fjet);
 250   fTree->SetBranchAddress("fjetaw", fjetaw);
 251   fTree->SetBranchAddress("fjphiw", fjphiw);
 252   fTree->SetBranchAddress("fjetal", fjetal);
 253   fTree->SetBranchAddress("fjphil", fjphil);
 254   // eT in EMCAL
 255   fTree->SetBranchAddress("fncell", &fncell);
 256   fTree->SetBranchAddress("fidcell", fidcell);
 257   fTree->SetBranchAddress("fetcell", fetcell);
 258   // eT in EMCAL grid for JF
 259   fTree->SetBranchAddress("fngrid", &fngrid);
 260   fTree->SetBranchAddress("fidgrid", fidgrid);
 261   fTree->SetBranchAddress("fetgrid", fetgrid);
 262   // 28-jan-2003
 263   fTree->SetBranchAddress("fnchp", &fnchp);
 264   fTree->SetBranchAddress("fpid", fpid);
 265   fTree->SetBranchAddress("fppt", fppt);
 266   fTree->SetBranchAddress("fpeta", fpeta);
 267   fTree->SetBranchAddress("fpphi", fpphi);
 268
 269   return kTRUE;
 270 }
 271
 272 void AliEMCALJetMicroDst::Print(Option_t* option) const
 273 {
 274         //
 275   if(option);
 276   if(fFile) {
 277     fFile->Print();
 278     if(fTree) fTree->Print();
 279     else Info("Print", "TRee is zero\n");
 280   } else {
 281      Info("Print", "File with TRee is closed \n Name of file %s(from fFileName", fFileName.Data());
 282   }
 283
 284   Info("Print", "******* Current(last) event *****");
 285   printf("#partons %2i \n", fnpart);
 286   for(Int_t i=0; i<fnpart; i++){
 287     printf("     %1i Pt %7.1f eta  %7.4f phi  %7.4f \n",
 288     i, fxpt[i],fxeta[i],fxphi[i]);
 289   }
 290   printf("#jets    %2i \n", fnjet);
 291   for(Int_t i=0; i<fnjet; i++){
 292     printf("     %1i Et %7.1f etaw %7.4f phiw %7.4f \n",
 293     i,fjet[i],fjetaw[i],fjphiw[i]);
 294   }
 295   printf(" Title %s \n", GetTitle());
 296 }
 297
 298 void AliEMCALJetMicroDst::Fill(AliRun *run, AliEMCALJetFinder* jetFinder, Int_t modeFilling)
 299 {
 300   //  modeFilling >=1 - fill info for EMCAL grid
 301   if(!run) run = gAlice;
 302   AliGenEventHeader* evHeader = run->GetHeader()->GenEventHeader();
 303   TString tmp(evHeader->ClassName());
 304   if(tmp.Contains("Hijing")) {
 305      AliGenHijingEventHeader *hijEvHeader = (AliGenHijingEventHeader*)evHeader;
 306      FillPartons(hijEvHeader);
 307   } else if(tmp.Contains("Pythia")) {
 308      FillPartons();
 309   } else {
 310     Error("Fill", "Wrong type of generator -> %s \n Info about partons will be absent",tmp.Data());
 311   }
 312
 313   FillJets(jetFinder);
 314
 315   if(modeFilling >= 1) {
 316      FillEtForEMCAL(jetFinder);
 317      FillEtForGrid(jetFinder);
 318      FillChargeParticles(jetFinder);
 319   } else {
 320      fncell = 0; // 27-jan-2003
 321      fngrid = 0; // 27-jan-2003
 322      fnchp  = 0;
 323      // negative - no signal
 324      fdecone = -1.;
 325      fptcone = -1.;
 326   }
 327
 328   FillJetsControl();  //24-jan-2003
 329
 330   fTree->Fill();
 331 }
 332
 333 void AliEMCALJetMicroDst::FillPartons(AliGenHijingEventHeader *header)
 334 {
 335   //Make partons arrays
 336   TLorentzVector parton[4];
 337   header->GetJets(parton[0], parton[1], parton[2], parton[3]);
 338
 339   fnpart = 4; //
 340   for(Int_t i=0; i<4; i++){
 341      fxpt[i]  = parton[i].Pt();
 342      fxeta[i] = parton[i].Eta();
 343      fxphi[i] = parton[i].Phi();
 344   }
 345 }
 346
 347 void AliEMCALJetMicroDst::FillPartons()
 348 {
 349   // for case of Pythia -> get info from full event record
 350
 351   fnpart = 2;
 352   TParticle *mPart;
 353   Int_t ind;
 354   for(Int_t i=6; i<8; i++){
 355      mPart = gAlice->GetMCApp()->Particle(i);
 356      ind   = i-6;
 357      fxpt[ind]  = mPart->Pt();
 358      fxeta[ind] = mPart->Eta();
 359      fxphi[ind] = mPart->Phi();
 360   }
 361 }
 362
 363 void AliEMCALJetMicroDst::FillJets(AliEMCALJetFinder* jetFinder)
 364 {
 365   // Fill Jets
 366   fnjet = 0;
 367   if(fDebug>1) Info("FillJets", "Debug");
 368   if(!jetFinder) {
 369     if(fDebug>1) Info("FillJets", "jetFinder is zero");
 370     return;
 371   }
 372   fnjet = jetFinder->Njets();
 373   if(fnjet>10) {
 374     if(fDebug>1) Warning("FillJets", "wrong value of jetFinder->Njets() %i ", fnjet);
 375     fnjet = 10;
 376   }
 377   //  fhNJet->Fill(njet);
 378   if(fDebug>1) Info("FillJets", "njet %i", fnjet);
 379   if(fnjet){
 380     for(Int_t i=0; i<fnjet; i++){
 381       fjet[i]   = jetFinder->JetEnergy(i);
 382       fjetaw[i] = jetFinder->JetEtaW(i);
 383       fjphiw[i] = jetFinder->JetPhiW(i);
 384       fjetal[i] = jetFinder->JetEtaL(i);
 385       fjphil[i] = jetFinder->JetPhiL(i);
 386     }
 387   }
 388 }
 389
 390 void AliEMCALJetMicroDst::FillEtForEMCAL(AliEMCALJetFinder* jetFinder)
 391 {
 392   // Fill Et for EMCAL
 393    fncell = 0;
 394    TH2F *hid = jetFinder->GetLegoEMCAL();
 395    if(!hid) return;
 396
 397    Double_t de = 0.;
 398    Int_t neta = hid->GetNbinsX(), nphi = hid->GetNbinsY();
 399    for(Int_t ieta=1; ieta<=neta; ieta++) {
 400       for(Int_t iphi=1; iphi<=nphi; iphi++) {
 401          de = hid->GetBinContent(ieta,iphi);
 402          if(de > 0.0) {
 403            fetcell[fncell] = Float_t(de);
 404            fidcell[fncell] = nphi*(ieta-1) + iphi;
 405            fncell++;
 406            if(fncell >= 13824) break;
 407            // Info("FillEtForEMCAL", " ncell %i6 id %i6 de %f \n", ncell, idcell[ncell], etcell[ncell]);
 408          }
 409       }
 410    }
 411    if(fnjet == 1) {
 412      // jet energy calculate around LP direction !!! - 10-mar-2003
 413       fdecone = jetFinder->EMCALConeEnergy(fjetal[0],fjphil[0]);
 414       fptcone = jetFinder->TrackConeEnergy(fjetal[0],fjphil[0]); // get from lego plot fo ch.part
 415       Info("FillEtForEMCAL", " njet %i Emcal in cone %f pt ch.part in cone %f\n", fnjet, fdecone, fptcone);
 416       Info("FillEtForEMCAL", " jet - decone - ptcone : %9.2f\n", fjet[0]-fdecone-fptcone);
 417    } else {
 418       fdecone = -1.;
 419       fptcone = -1.;
 420    }
 421
 422    Info("FillEtForEMCAL", "neta %3i nphi %3i # array size %i Sum.Et %f\n",
 423    neta,nphi, fncell, hid->Integral());
 424 }
 425
 426 void AliEMCALJetMicroDst::FillEtForGrid(AliEMCALJetFinder* jetFinder)
 427 {
 428   // Fill ET for Grid
 429    TH2F *hid = jetFinder->GetLego();
 430    if(!hid) return;
 431
 432    FillArrays(hid, fngrid, fidgrid, fetgrid);
 433 }
 434
 435 void AliEMCALJetMicroDst::FillArrays(TH2* hid, Int_t &n, Int_t *id, Float_t *et)
 436 {
 437   // Fill arays
 438    n = 0;
 439    Double_t de = 0.;
 440    Int_t neta = hid->GetNbinsX(), nphi = hid->GetNbinsY();
 441    for(Int_t ieta=1; ieta<=neta; ieta++) {
 442       for(Int_t iphi=1; iphi<=nphi; iphi++) {
 443          de = hid->GetBinContent(ieta,iphi);
 444          if(de > 0.0) {
 445            et[n] = Float_t(de);
 446            id[n] = nphi*(ieta-1) + iphi;
 447            n++;
 448            if(n >= 13824) break;
 449          }
 450       }
 451    }
 452    Info("FillArrays", "neta %3i nphi %3i # array size %i Sum.Et %f\n",
 453    neta, nphi, n, hid->Integral());
 454 }
 455
 456 void AliEMCALJetMicroDst::FillChargeParticles(AliEMCALJetFinder* jetFinder)
 457 {
 458   // 28-jan-2003 for fullness ; 18-mar - sometimes
 459   fnchp = 0;
 460   Int_t gid=0;
 461   for(Int_t i=0; i<jetFinder->fNt; i++) {
 462     //     fPdgT[i];
 463      if(jetFinder->fTrackList[i] >= 1) {
 464         sgpdge_(jetFinder->fPdgT[i], gid);
 465         fpid[fnchp]  = gid;
 466         fppt[fnchp]  = jetFinder->fPtT[i];
 467         fpeta[fnchp] = jetFinder->fEtaT[i];
 468         fpphi[fnchp] = jetFinder->fPhiT[i];
 469         fnchp++;
 470      }
 471   }
 472   Info("FillChargedParticles", "fNtS %i : nchp %i -> %i\n", jetFinder->fNtS, fnchp, jetFinder->fNtS - fnchp);
 473 }
 474
 475 void AliEMCALJetMicroDst::FillJetsControl()
 476 {
 477   // see FillJets(AliEMCALJetFinder* jetFinder) and FillPartons
 478   fhNJet->Fill(fnjet);
 479   for(Int_t i=0; i<fnjet; i++){
 480      fhPtJet->Fill(fjet[i]);
 481      fhEtaPhiJet->Fill(fjetaw[i],fjphiw[i]);
 482   }
 483
 484   for(Int_t i=0; i < fnpart; i++){
 485      fhEtaPhiPart->Fill(fxeta[i], fxphi[i]);
 486      fhPtPart->Fill(fxpt[i]);
 487   }
 488
 489   Double_t sum = 0.0;
 490   fhNcell->Fill(fncell);
 491   for(Int_t i=0; i < fncell; i++){
 492     fhCellId->Fill(fidcell[i]);
 493     fhCellEt->Fill(fetcell[i]);
 494     sum += Double_t(fetcell[i]);
 495   }
 496   fhSumEt->Fill(sum);
 497
 498   sum = 0.0;
 499   fhNgrid->Fill(fngrid);
 500   for(Int_t i=0; i < fngrid; i++){
 501     fhGridId->Fill(fidgrid[i]);
 502     fhGridEt->Fill(fetgrid[i]);
 503     sum += Double_t(fetgrid[i]);
 504   }
 505   fhSumEtGrForJF->Fill(sum);
 506 }
 507
 508 Int_t AliEMCALJetMicroDst::GetEntry(Int_t entry)
 509 {
 510   // Read contents of entry.
 511    if (!fTree) {
 512       Error("GetEntry", "define TTree");
 513       return -1;
 514    }
 515    return fTree->GetEntry(entry);
 516 }
 517
 518 Bool_t AliEMCALJetMicroDst::GetParton(Int_t i, Float_t& pt, Float_t& eta, Float_t& phi) const
 519 {
 520   // Get parton
 521   if(i>=0 && i<fnpart) {
 522     pt  = fxpt[i];
 523     eta = fxeta[i];
 524     phi = fxphi[i];
 525     return kTRUE;
 526   } else return kFALSE;
 527 }
 528
 529 Bool_t AliEMCALJetMicroDst::GetParton(Int_t i, TVector3& vec) const
 530 {
 531   // Get Parton
 532   static Float_t pt, eta, phi;
 533
 534   if(!GetParton(i, pt, eta, phi)) return kFALSE;
 535
 536   FillVector(pt, eta, phi, vec);
 537   return kTRUE;
 538 }
 539
 540 Bool_t AliEMCALJetMicroDst::GetJet(Int_t i,Int_t mode,Float_t& pt, Float_t& eta, Float_t& phi) const
 541 {
 542   // mode=1(W) mode=any(L)
 543   if(i>=0 && i<fnjet) {
 544     pt  = fjet[i];
 545     if(mode==1) {
 546       eta = fjetaw[i];
 547       phi = fjphiw[i];
 548     } else {
 549       eta = fjetal[i];
 550       phi = fjphil[i];
 551     }
 552     return kTRUE;
 553   } else return kFALSE;
 554 }
 555
 556 Bool_t AliEMCALJetMicroDst::GetJet(Int_t i, Int_t mode, TVector3& vec) const
 557 {
 558   // Get Jet
 559   static Float_t pt, eta, phi;
 560
 561   if(!GetJet(i, mode, pt, eta, phi)) return kFALSE;
 562   FillVector(pt, eta, phi, vec);
 563   return kTRUE;
 564 }
 565
 566 void AliEMCALJetMicroDst::Test()
 567 {
 568   // Test
 569   if(!fFile || !fTree ) {
 570     Info("Test", "define file with proper TTree !");
 571     return;
 572   }
 573   Int_t nbytes=0, nb=0, nentries=Int_t(fTree->GetEntries());
 574   for(Int_t i=0; i<nentries; i++){
 575     nb = fTree->GetEntry(i);
 576     nbytes += nb;
 577     for(Int_t j=0; j<fnpart; j++){
 578       fhEtaPhiPart->Fill(fxeta[j], fxphi[j]);
 579       fhPtPart->Fill(fxpt[j]);
 580     }
 581
 582     fhNJet->Fill(fnjet);
 583     if(fnjet){
 584       for(Int_t j=0; j<fnjet; j++) {
 585         fhPtJet->Fill(fjet[j]);
 586         fhEtaPhiJet->Fill(fjetaw[j],fjphiw[j]);
 587       }
 588     }
 589   }
 590   Info("Test", "Entries %5i Bytes %10i\n", nentries, nbytes);
 591 }
 592
 593 void AliEMCALJetMicroDst::FillVector(Float_t pt, Float_t eta, Float_t phi, TVector3& vec)
 594 {
 595   // service function
 596   static Float_t px, py, pz;
 597
 598   px = pt*TMath::Cos(phi);
 599   py = pt*TMath::Sin(phi);
 600   pz = pt*TMath::SinH(eta);  // sinh(eta) = cot(theta)
 601
 602   vec.SetXYZ(px, py, pz);
 603 }
 604
 605 void AliEMCALJetMicroDst::GetEtaPhi(Int_t id, Double_t &eta, Double_t &phi) const
 606 {
 607   // see AliEMCALGeometry
 608   static Int_t ieta, iphi, nphi=144, neta=96;
 609   static Double_t phiMax=(2.0/3.0)*TMath::Pi(), phiMin=0.0;
 610   static Double_t phiStep=(phiMax-phiMin)/nphi, phiBeg = phiMin + phiStep/2.;
 611   static Double_t etaMax=0.7, etaMin=-etaMax;
 612   static Double_t etaStep=(etaMax-etaMin)/neta, etaBeg = etaMin + etaStep/2.;
 613
 614   ieta = (id-1)/nphi + 1; // id = nphi*(ieta-1) + iphi
 615   iphi = id - nphi*(ieta-1);
 616   if(ieta<=0 || ieta>neta) {
 617     Fatal("GetEtaPhi", "wrong id %i(ieta %i,iphi %i) : nphi %i neta %i", id,iphi,ieta, nphi,neta);
 618
 619   }
 620
 621   eta  = etaBeg + etaStep*(ieta-1);
 622   phi  = phiBeg + phiStep*(iphi-1);
 623 }
 624
 625 TVector3& AliEMCALJetMicroDst::GetCellVector(Int_t i) const
 626 {
 627   // Get cell vector
 628   static Double_t eta,phi;
 629   static TVector3 vec;
 630   vec.SetXYZ(0.,0.,0.);
 631   if(i>=0 && i<fncell) {
 632      GetEtaPhi(fidcell[i], eta, phi);
 633      vec.SetPtEtaPhi(Double_t(fetcell[i]),eta,phi);
 634   }
 635   return vec;
 636 }
 637
 638 TVector3& AliEMCALJetMicroDst::GetGridVector(Int_t i) const
 639 {
 640   // Get grid vector
 641   static Double_t eta,phi;
 642   static TVector3 vec;
 643   vec.SetXYZ(0.,0.,0.);
 644   if(i>=0 && i<fngrid) {
 645      GetEtaPhi(fidgrid[i], eta, phi);
 646      vec.SetPtEtaPhi(Double_t(fetgrid[i]),eta,phi);
 647   }
 648   return vec;
 649 }
 650
 651 Double_t AliEMCALJetMicroDst::GetSumInCone(TVector3 &jet,Int_t nc, Float_t *et,Float_t *eta,Float_t *phi, Double_t cellEtCut, Double_t rJet) const
 652 {
 653   // Get Sum in cone
 654   static Double_t sum=0.;
 655   static TVector3 cell(0., 0., 0.);
 656   if(nc<=0 || et==0 || eta==0 || phi==0) {
 657     Error("GetSumInCone", "nc %d %f %f %f ", nc, et, eta, phi);
 658     return -1.;
 659   }
 660
 661   sum=0.;
 662   //  jet.SetPtEtaPhi(jet[0],jetaw[0],jphiw[0]); // must be one jet !!
 663   Info("GetSumInCone", "jet.Mag() %f : njet %i\n", jet.Mag(), fnjet);
 664   for(Int_t i=0; i<nc; i++){
 665     if(et[i] < cellEtCut)    continue;
 666     cell.SetPtEtaPhi(et[i], eta[i], phi[i]);
 667     if(jet.DeltaR(cell) > rJet) continue;
 668     sum += et[i];
 669   }
 670   Info("GetSumCone", "Sum %f \n", sum);
 671   return sum;
 672 }
 673
 674 Double_t AliEMCALJetMicroDst::GetEmcalEtInCone(TVector3 &jet, Double_t cellEtCut, Double_t rJet)
 675 {
 676   // Get EMCAL Et in cone
 677   Int_t nc = fncell;
 678   if(nc<=0) return 0.;
 679
 680   Float_t *et=fetcell, *eta=new Float_t[nc], *phi=new Float_t[nc];
 681   Double_t etaCell=0., phiCell=0., eTotal=0;
 682
 683   for(Int_t i=0; i<nc; i++) {
 684     GetEtaPhi(fidcell[i], etaCell, phiCell);
 685     eta[i] = etaCell;
 686     phi[i] = phiCell;
 687   }
 688
 689   eTotal = GetSumInCone(jet, nc, et,eta,phi, cellEtCut,rJet);
 690   delete [] eta;
 691   delete [] phi;
 692
 693   return eTotal;
 694 }
 695
 696 Double_t AliEMCALJetMicroDst::GetTpcPtInCone(TVector3 &jet,Double_t cellEtCut, Double_t rJet)
 697 {
 698   // Get TPC PT in cone
 699   if(fnchp<=0) return 0.;
 700   return GetSumInCone(jet, fnchp, fppt,fpeta,fpphi, cellEtCut,rJet);
 701 }
 702
 703 Double_t AliEMCALJetMicroDst::GetSum(Int_t n, Float_t *ar, Double_t cut) const
 704 {
 705   // 25-apr-2003
 706   Double_t sum=0.0;
 707   if(n<=0 || ar==0) return sum;
 708   for(Int_t i=0; i<n; i++) {if(ar[i]>=cut) sum += Double_t(ar[i]);}
 709   return sum;
 710 }
 711
 712 void AliEMCALJetMicroDst::Close()
 713 {
 714   // Close
 715   fFile->Write();
 716   fTree->Print();
 717   fFile->Close();
 718   fFile = 0;
 719   fTree = 0;
 720 }
 721
 722 void AliEMCALJetMicroDst::Browse(TBrowser* b) const
 723 {
 724   // Browse
 725    if(fTree)      b->Add((TObject*)fTree);
 726    if(fListHist)  b->Add((TObject*)fListHist);
 727    //   TObject::Browse(b);
 728 }
 729
 730 Bool_t  AliEMCALJetMicroDst::IsPythiaDst() const
 731 {
 732   // Is Pythia DST
 733   TString st(GetTitle());
 734   if(st.Contains("py",TString::kIgnoreCase)||!st.Contains("hijing", TString::kIgnoreCase)) return kTRUE;
 735   else return kFALSE;
 736 }
 737
 738 Bool_t  AliEMCALJetMicroDst::IsFolder() const
 739 {
 740   // Is folder
 741   if(fTree || fListHist) return kTRUE;
 742   else                   return kFALSE;
 743 }
 744
 745 TList* AliEMCALJetMicroDst::MoveHistsToList(const char* name, Bool_t putToBrowser)
 746 {
 747   // Move HIST to list
 748   gROOT->cd();
 749   TIter nextHist(gDirectory->GetList());
 750   TList *list = new TList;
 751   list->SetName(name);
 752   TObject *objHist;
 753   while((objHist=nextHist())){
 754     if (!objHist->InheritsFrom("TH1")) continue;
 755     ((TH1*)objHist)->SetDirectory(0); // Remove from gROOT
 756     list->Add(objHist);
 757   }
 758   if(putToBrowser) gROOT->GetListOfBrowsables()->Add((TObject*)list);
 759   return list;
 760 }