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