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