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