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