]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Last version
authorpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Jul 2003 16:24:19 +0000 (16:24 +0000)
committerpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Jul 2003 16:24:19 +0000 (16:24 +0000)
EMCAL/AliEMCALJetMicroDst.cxx
EMCAL/AliEMCALJetMicroDst.h

index c79814b425c29646ad3c066a302513598a144f51..35bcb83f3eb288365fec8136fe5e47f7f291bc8c 100644 (file)
@@ -20,6 +20,8 @@
 //*
 
 #include "AliEMCALJetMicroDst.h"
+#include <assert.h>
+#include <iostream>
 #include "AliRun.h"
 #include "AliHeader.h"
 #include "AliGenEventHeader.h"
@@ -37,25 +39,40 @@ TString nameTree("jetMDST"); // 7-feb-2002
 
 TH1F*  hPtPart, *hNJet, *hPtJet;
 TH2F*  hEtaPhiPart, *hEtaPhiJet;
+TH1F*  hNcell, *hCellId, *hCellEt, *hSumEt;
+TH1F*  hNgrid, *hGridId, *hGridEt, *hSumEtGrForJF;
+
+extern "C" void sgpdge_(Int_t &i, Int_t &pdggea);
 
 AliEMCALJetMicroDst::AliEMCALJetMicroDst(char *name, char *tit) : TNamed(name,tit)
 {
   fFile  = 0;
   fTree  = 0;
   fDebug = 0;
-  fName  = "dstTest.root"; // default name
     
 //  Don't add histos to the current directory
 //  TH1::AddDirectory(0);
   gROOT->cd();
   hPtPart     = new TH1F("hPtPart","P_{T} for partons", 300, 0., 300.);
+  // 16-jan-2002 - new limit fo phi 
   hEtaPhiPart = new TH2F("hEtaPhiPart","#eta #phi distr.for partons after HSc", 
-                      28, -0.7, 0.7, 21, 0.0, TMath::Pi()*2./3.);
+                      28, -0.7, 0.7, 21, TMath::Pi()/3., TMath::Pi());
 
   hNJet  = new TH1F("hNJet","number of jets", 11, -0.5, 10.5);
-  hPtJet = new TH1F("hPtJet","P_{T} for jets", 300, 0., 300.);
+  hPtJet = new TH1F("hPtJet","P_{T} for jets", 500, 0., 500.);
   hEtaPhiJet = new TH2F("hEtaPhiJet","#eta #phi distr.for jets (W)", 
-                      28, -0.7, 0.7, 21, 0.0, TMath::Pi()*2./3.);
+                      28, -0.7, 0.7, 21, TMath::Pi()/3., TMath::Pi());
+
+  hNcell  = new TH1F("hNcell","#cell with de>0.0 for EMCAL", 1400, 0.0, 14000.);
+  hCellId = new TH1F("hCellId","cell ID with de>0.0 for EMCAL", 1400, 0.0, 14000.);
+  hCellEt = new TH1F("hCellEt","cell Et for EMCAL", 1000, 0.0, 10.);
+  hSumEt  = new TH1F("hSumEt","sum Et for EMCAL", 1000, 0.0, 1000.);
+
+  hNgrid  = new TH1F("hNgrid","#cell with de>0.0 in EMCAL grid for JF", 1400, 0.0, 14000.);
+  hGridId = new TH1F("hGridId","cell ID with de>0.0 in EMCAL grid for JF", 1400, 0.0, 14000.);
+  hGridEt = new TH1F("hGridEt","cell Et in EMCAL grid for JF", 1000, 0.0, 10.);
+  hSumEtGrForJF  = new TH1F("hSumEtGrForJF","sum Et in EMCAL grid for JF", 1000, 0.0, 1000.);
+
   fListHist = MoveHistsToList("Hist For AliEMCALJetMicroDst", kFALSE); 
 }
 
@@ -70,21 +87,43 @@ Bool_t AliEMCALJetMicroDst::Create(TFile *file)
     printf("<E> AliEMCALJetMicroDst::create -> define TFile for output\n");
     return kFALSE;
   }
-  fFile = file;
+  fFile     = file;
+  fFileName = fFile->GetName();
   fFile->cd();
   fTree = new TTree(nameTree.Data(),"Temporary micro DST for jets analysis");
+  // for jet calibration - 4-mar-2003
+  fTree->Branch("decone", &decone, "decone/F");
+  fTree->Branch("ptcone", &ptcone, "ptcone/F");
   // partons 
-  fTree->Branch("npart", &fNpart, "npart/I");
-  fTree->Branch("xpt",  fXpt,  "xpt[npart]/F");
-  fTree->Branch("xeta", fXeta, "xeta[npart]/F");
-  fTree->Branch("xphi", fXphi, "xphi[npart]/F");
+  fTree->Branch("npart", &npart, "npart/I");
+  fTree->Branch("xpt",  xpt,  "xpt[npart]/F");
+  fTree->Branch("xeta", xeta, "xeta[npart]/F");
+  fTree->Branch("xphi", xphi, "xphi[npart]/F");
   // jets 
-  fTree->Branch("njet", &fNjet, "njet/I");
-  fTree->Branch("jet",  fJet,  "jet[njet]/F");
-  fTree->Branch("jetaw", fJetaW, "jetaw[njet]/F");
-  fTree->Branch("jphiw", fJphiW, "jphiw[njet]/F");
-  fTree->Branch("jetal", fJetaL, "jetal[njet]/F");
-  fTree->Branch("jphil", fJphiL, "jphil[njet]/F");
+  fTree->Branch("njet", &njet, "njet/I");
+  fTree->Branch("jet",  jet,  "jet[njet]/F");
+  fTree->Branch("jetaw", jetaw, "jetaw[njet]/F");
+  fTree->Branch("jphiw", jphiw, "jphiw[njet]/F");
+  fTree->Branch("jetal", jetal, "jetal[njet]/F");
+  fTree->Branch("jphil", jphil, "jphil[njet]/F");
+
+  // Et in EMCAL itself 
+  fTree->Branch("ncell", &ncell, "ncell/I");
+  fTree->Branch("idcell", idcell, "idcell[ncell]/I");
+  fTree->Branch("etcell", etcell, "etcell[ncell]/F");
+
+  // Et in EMCAL grid for JF 
+  fTree->Branch("ngrid", &ngrid,  "ngrid/I");
+  fTree->Branch("idgrid", idgrid, "idgrid[ngrid]/I");
+  fTree->Branch("etgrid", etgrid, "etgrid[ngrid]/F");
+
+  // charge particle which hit to EMCAL
+  fTree->Branch("nchp", &nchp, "nchp/I");
+  fTree->Branch("pid", pid, "pid[nchp]/I");
+  fTree->Branch("ppt", ppt, "ppt[nchp]/F");
+  fTree->Branch("peta", peta, "peta[nchp]/F");
+  fTree->Branch("pphi", pphi, "pphi[nchp]/F");
+
   return kTRUE;
 }
 
@@ -92,9 +131,9 @@ Bool_t AliEMCALJetMicroDst::Create(const char *fname)
 {
   TFile *file = new TFile(fname, "RECREATE");
   if(file) {
-    fName = fname;
+    //    fNameFile = fname;
     return Create(file);
-  } else     return kFALSE;
+  } else   return kFALSE;
 }
 
 Bool_t AliEMCALJetMicroDst::Open(const char *fname)
@@ -103,32 +142,120 @@ Bool_t AliEMCALJetMicroDst::Open(const char *fname)
   if(strlen(fname)) fName = fname;
   TFile *file = new TFile(fName.Data(), "READ");
   if(file) {
-    printf("<I> open file %s \n",fName.Data()); 
-    return Initialize(file);
+    Bool_t ini =  Initialize(file);
+    printf("<I> open file %s : initialize TTree %i \n",fName.Data(), Int_t(ini)); 
+    return ini;
   } else {
     printf("<E> can not open file %s \n",fName.Data()); 
     return kFALSE;
   }
 }
 
+const Char_t* AliEMCALJetMicroDst::DefineName(const Int_t mode)
+{
+  static TString dir, name;
+  //  dir = "jetDST/"; // 24-jan-2003
+  dir = "/auto/alice/pavlinov/jet/microDST/"; // 24-jan-2003
+  switch (mode) {
+  case 1: // for characteristic of BG
+    name = dir + "mDst1_1.root"; // Bg 2000 - first version 
+    SetTitle("Bg2000");
+    break;
+  case 2: // for characteristic of BG
+    name = dir + "mDst2_1.root"; // Bg 4000 - first version 
+    SetTitle("Bg4000");
+    break;
+  case 3: // for characteristic of BG
+    name = dir + "mDst3_1.root"; // Bg 8000 - first version 
+    SetTitle("Bg8000");
+    break;
+  case 4: // Central Hijing - 18-mar-2003
+    name  = dir  + "march/";
+    name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0HijingCentral.root";
+    SetTitle("HijingCentral");
+    break;
+  case 5: // Para Hijing (Dn/Dy=8000) - 21-mar-2003
+    name  = dir  + "march/";
+    name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0ParaHijing8000.root";
+    SetTitle("HIJINGparaDnDy8000");
+    break;
+  case 6: // Para Hijing (Dn/Dy=4000) - 21-mar-2003
+    name  = dir  + "march/";
+    name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0ParaHijing4000.root";
+    SetTitle("HIJINGparaDnDy4000");
+    break;
+  case 7: // Para Hijing (Dn/Dy=2000) - 21-mar-2003
+    name  = dir  + "march/";
+    name += "jF_R0.50MinCell1.0PtCut0.0EtSeed4.0MinEt40.0BGSubtr0SF11.6Smear0Eff0ParaHijing2000.root";
+    SetTitle("HIJINGparaDnDy2000");
+    break;
+  case 11: // pure PYTHIA with default value of parameters
+    name = dir + "jF_R0.50MinCell0.0PtCut0.0EtSeed8.0MinEt40.0BGSubtr0SF11.6.root";
+    break;
+  case 12: // 0 + background 
+    name = dir + "jF_R0.50MinCell0.0PtCut0.0EtSeed8.0MinEt40.0BGSubtr0SF11.6kBackground2000.root";
+    break;
+    // calibration case 
+  case 101:
+    name  = dir  + "march/";
+    name += "Pythia100_1.root";
+    SetTitle("Pythia100_1");
+    break;
+  case 102: // 2-apr-2003
+    name  = dir  + "march/";
+    name += "Pythia50_1.root";
+    SetTitle("Pythia50_1");
+    //    name = "microDst3th.root"; // 101 + (smearing and eff) - 14-mar-2003
+    break;
+  case 103:// 4-apr-2003
+    name  = dir  + "march/";
+    name += "Pythia200_1.root";
+    SetTitle("Pythia200_1");
+    //    name = "microDst4th.root"; // 102 + MinCell= 1.0
+    break;
+  default:
+    printf("AliEMCALJetMicroDst::DefineName : NO  D E F A U L T : mode %i\n", mode);
+    assert(0);
+  }
+  printf("mode %5i file : %s : Title %s\n", mode, name.Data(), GetTitle());
+  return name.Data();
+}
+
 Bool_t AliEMCALJetMicroDst::Initialize(TFile *file)
 {
   if(file) fFile = file; 
   fFile->cd();
   fTree = (TTree*)fFile->Get(nameTree.Data());
   if(!fTree) return kFALSE;
+  // for jet calibration - 4-mar-2003
+  fTree->SetBranchAddress("decone",&decone);
+  fTree->SetBranchAddress("ptcone",&ptcone);
   // partons
-  fTree->SetBranchAddress("npart",&fNpart);
-  fTree->SetBranchAddress("xpt",   fXpt);
-  fTree->SetBranchAddress("xeta",  fXeta);
-  fTree->SetBranchAddress("xphi",  fXphi);
+  fTree->SetBranchAddress("npart",&npart);
+  fTree->SetBranchAddress("xpt",   xpt);
+  fTree->SetBranchAddress("xeta",  xeta);
+  fTree->SetBranchAddress("xphi",  xphi);
   // jets 
-  fTree->SetBranchAddress("njet", &fNjet);
-  fTree->SetBranchAddress("jet",   fJet);
-  fTree->SetBranchAddress("jetaw", fJetaW);
-  fTree->SetBranchAddress("jphiw", fJphiW);
-  fTree->SetBranchAddress("jetal", fJetaL);
-  fTree->SetBranchAddress("jphil", fJphiL);
+  fTree->SetBranchAddress("njet", &njet);
+  fTree->SetBranchAddress("jet",   jet);
+  fTree->SetBranchAddress("jetaw", jetaw);
+  fTree->SetBranchAddress("jphiw", jphiw);
+  fTree->SetBranchAddress("jetal", jetal);
+  fTree->SetBranchAddress("jphil", jphil);
+  // eT in EMCAL
+  fTree->SetBranchAddress("ncell", &ncell);
+  fTree->SetBranchAddress("idcell", idcell);
+  fTree->SetBranchAddress("etcell", etcell);
+  // eT in EMCAL grid for JF
+  fTree->SetBranchAddress("ngrid", &ngrid);
+  fTree->SetBranchAddress("idgrid", idgrid);
+  fTree->SetBranchAddress("etgrid", etgrid);
+  // 28-jan-2003
+  fTree->SetBranchAddress("nchp", &nchp);
+  fTree->SetBranchAddress("pid", pid);
+  fTree->SetBranchAddress("ppt", ppt);
+  fTree->SetBranchAddress("peta", peta);
+  fTree->SetBranchAddress("pphi", pphi);
 
   return kTRUE;
 }
@@ -140,24 +267,27 @@ void AliEMCALJetMicroDst::Print(Option_t* option) const
     fFile->Print();
     if(fTree) fTree->Print();
     else printf("<I> TRee is zero\n");
-  } else printf("<I> File with TRee is zero\n");
-  printf("Default name of file %s\n", fName.Data());
-
-  printf("\n #partons %2i ", fNpart);
-  for(Int_t i=0; i<fNpart; i++){
-    printf("\n     %1i Pt %7.1f eta  %f7.4 phi  %f7.4 ",
-    i, fXpt[i],fXeta[i],fXphi[i]);  
+  } else {
+     printf("<I> File with TRee is closed \n");
+     printf(" Name of file %s(from fFileName \n", fFileName.Data());
   }
 
-  printf("\n #jets    %2i ", fNjet);
-  for(Int_t i=0; i<fNjet; i++){
-    printf("\n     %1i Et %7.1f etaW %f7.4 phiW %f7.4 ",
-    i,fJet[i],fJetaW[i],fJphiW[i]);  
+  printf("******* Current(last) event ***** \n");
+  printf("#partons %2i \n", npart);
+  for(Int_t i=0; i<npart; i++){
+    printf("     %1i Pt %7.1f eta  %7.4f phi  %7.4f \n",
+    i, xpt[i],xeta[i],xphi[i]);  
   }
+  printf("#jets    %2i \n", njet);
+  for(Int_t i=0; i<njet; i++){
+    printf("     %1i Et %7.1f etaw %7.4f phiw %7.4f \n",
+    i,jet[i],jetaw[i],jphiw[i]);  
+  }
+  printf(" Title %s \n", GetTitle());
 }
 
-void AliEMCALJetMicroDst::Fill(AliRun *run, AliEMCALJetFinder* jetFinder)
-{
+void AliEMCALJetMicroDst::Fill(AliRun *run, AliEMCALJetFinder* jetFinder, Int_t modeFilling)
+{//  modeFilling >=1 - fill info for EMCAL grid
   if(!run) run = gAlice;
   AliGenEventHeader* evHeader = run->GetHeader()->GenEventHeader();
   TString tmp(evHeader->ClassName());
@@ -167,11 +297,27 @@ void AliEMCALJetMicroDst::Fill(AliRun *run, AliEMCALJetFinder* jetFinder)
   } else if(tmp.Contains("Pythia")) {
      FillPartons();     
   } else {
-    printf("\n <E> Wrong type of generator -> %s ",tmp.Data()); 
-    printf("\n     Inof about partons will be absent "); 
+    printf("<E> Wrong type of generator -> %s \n",tmp.Data()); 
+    printf("     Info about partons will be absent \n"); 
   }
 
   FillJets(jetFinder);
+
+  if(modeFilling >= 1) {
+     FillEtForEMCAL(jetFinder);
+     FillEtForGrid(jetFinder);
+     FillChargeParticles(jetFinder);
+  } else {
+     ncell = 0; // 27-jan-2003
+     ngrid = 0; // 27-jan-2003
+     nchp  = 0;
+     // negative - no signal 
+     decone = -1.;
+     ptcone = -1.;
+  }
+
+  FillJetsControl();  //24-jan-2003
+
   fTree->Fill();
 }
 
@@ -180,64 +326,168 @@ void AliEMCALJetMicroDst::FillPartons(AliGenHijingEventHeader *header)
   TLorentzVector parton[4];
   header->GetJets(parton[0], parton[1], parton[2], parton[3]);
 
-  fNpart = 4; // 
+  npart = 4; // 
   for(Int_t i=0; i<4; i++){
-     fXpt[i]  = parton[i].Pt();
-     fXeta[i] = parton[i].Eta();
-     fXphi[i] = parton[i].Phi();
-     if(i>1) {
-       hEtaPhiPart->Fill(parton[i].Eta(), parton[i].Phi());
-       hPtPart->Fill(parton[i].Pt());
-     }
+     xpt[i]  = parton[i].Pt();
+     xeta[i] = parton[i].Eta();
+     xphi[i] = parton[i].Phi();
   }
 }
 
 void AliEMCALJetMicroDst::FillPartons()
 {// for case of Pythia -> get info from full event record
 
-  fNpart = 2;
+  npart = 2;
   TParticle *MPart;
   Int_t ind;
   for(Int_t i=6; i<8; i++){
      MPart = gAlice->Particle(i);
      ind   = i-6;
-     fXpt[ind]  = MPart->Pt();
-     fXeta[ind] = MPart->Eta();
-     fXphi[ind] = MPart->Phi();
-
-     hEtaPhiPart->Fill(fXeta[ind], fXphi[ind]);
-     hPtPart->Fill(fXpt[ind]);
+     xpt[ind]  = MPart->Pt();
+     xeta[ind] = MPart->Eta();
+     xphi[ind] = MPart->Phi();
   }
 }
 
 void AliEMCALJetMicroDst::FillJets(AliEMCALJetFinder* jetFinder)
 {
-  fNjet = 0;
+  njet = 0;
   if(fDebug>1) printf("\n<I> AliEMCALJetMicroDst::FillJets"); 
   if(!jetFinder) {
     if(fDebug>1) printf("\n : jetFinder is zero"); 
     return;
   }
-  fNjet = jetFinder->Njets();
-  if(fNjet>10) {
-    if(fDebug>1) printf("\n <W> wrong value of jetFinder->Njets() %i ", fNjet); 
-    fNjet = 10;
+  njet = jetFinder->Njets();
+  if(njet>10) {
+    if(fDebug>1) printf("\n <W> wrong value of jetFinder->Njets() %i ", njet); 
+    njet = 10;
   }
-  hNJet->Fill(fNjet);
-  if(fDebug>1) printf("\n <I> fNjet %i", fNjet); 
-  if(fNjet){
-    for(Int_t i=0; i<fNjet; i++){
-      fJet[i]   = jetFinder->JetEnergy(i);
-      fJetaW[i] = jetFinder->JetEtaW(i);
-      fJphiW[i] = jetFinder->JetPhiW(i);
-      fJetaL[i] = jetFinder->JetEtaL(i);
-      fJphiL[i] = jetFinder->JetPhiL(i);
-      hPtJet->Fill(fJet[i]);
-      hEtaPhiJet->Fill(fJetaW[i],fJphiW[i]); 
+  //  hNJet->Fill(njet);
+  if(fDebug>1) printf("\n <I> njet %i", njet); 
+  if(njet){
+    for(Int_t i=0; i<njet; i++){
+      jet[i]   = jetFinder->JetEnergy(i);
+      jetaw[i] = jetFinder->JetEtaW(i);
+      jphiw[i] = jetFinder->JetPhiW(i);
+      jetal[i] = jetFinder->JetEtaL(i);
+      jphil[i] = jetFinder->JetPhiL(i);
     }
   }
 }
 
+void AliEMCALJetMicroDst::FillEtForEMCAL(AliEMCALJetFinder* jetFinder)
+{
+   ncell = 0;
+   TH2F *hid = jetFinder->GetLegoEMCAL();
+   if(!hid) return;
+
+   Double_t de = 0.;
+   Int_t neta = hid->GetNbinsX(), nphi = hid->GetNbinsY();
+   for(Int_t ieta=1; ieta<=neta; ieta++) {
+      for(Int_t iphi=1; iphi<=nphi; iphi++) {
+        de = hid->GetBinContent(ieta,iphi);
+         if(de > 0.0) {
+          etcell[ncell] = Float_t(de);
+           idcell[ncell] = nphi*(ieta-1) + iphi;
+          ncell++;
+           if(ncell >= 13824) break; 
+          //           printf(" ncell %i6 id %i6 de %f \n", ncell, idcell[ncell], etcell[ncell]); 
+         }
+      }
+   }
+   if(njet == 1) {
+     // jet energy calculate around LP direction !!! - 10-mar-2003 
+      decone = jetFinder->EMCALConeEnergy(jetal[0],jphil[0]); 
+      ptcone = jetFinder->TrackConeEnergy(jetal[0],jphil[0]); // get from lego plot fo ch.part
+      printf(" njet %i Emcal in cone %f pt ch.part in cone %f\n", njet, decone, ptcone); 
+      printf(" jet - decone - ptcone : %9.2f\n", jet[0]-decone-ptcone);
+   } else {
+      decone = -1.;
+      ptcone = -1.;
+   }
+
+   printf(" FillEtForEMCAL : neta %3i nphi %3i # array size %i Sum.Et %f\n", 
+   neta,nphi, ncell, hid->Integral());
+}
+
+void AliEMCALJetMicroDst::FillEtForGrid(AliEMCALJetFinder* jetFinder)
+{
+   TH2F *hid = jetFinder->GetLego();
+   if(!hid) return;
+
+   FillArrays(hid, ngrid, idgrid, etgrid);
+}
+
+void AliEMCALJetMicroDst::FillArrays(TH2* hid, Int_t &n, Int_t *id, Float_t *et)
+{
+   n = 0;
+   Double_t de = 0.;
+   Int_t neta = hid->GetNbinsX(), nphi = hid->GetNbinsY();
+   for(Int_t ieta=1; ieta<=neta; ieta++) {
+      for(Int_t iphi=1; iphi<=nphi; iphi++) {
+        de = hid->GetBinContent(ieta,iphi);
+         if(de > 0.0) {
+          et[n] = Float_t(de);
+           id[n] = nphi*(ieta-1) + iphi;
+          n++;
+           if(n >= 13824) break; 
+         }
+      }
+   }
+   printf(" AliEMCALJetMicroDst::FillArrays : neta %3i nphi %3i # array size %i Sum.Et %f\n", 
+   neta, nphi, n, hid->Integral());
+}
+
+void AliEMCALJetMicroDst::FillChargeParticles(AliEMCALJetFinder* jetFinder)
+{// 28-jan-2003 for fullness ; 18-mar - sometimes 
+  nchp = 0;
+  Int_t gid=0;
+  for(Int_t i=0; i<jetFinder->fNt; i++) {
+    //     fPdgT[i];
+     if(jetFinder->fTrackList[i] >= 1) {
+        sgpdge_(jetFinder->fPdgT[i], gid);
+        pid[nchp]  = gid; 
+        ppt[nchp]  = jetFinder->fPtT[i];
+        peta[nchp] = jetFinder->fEtaT[i];
+        pphi[nchp] = jetFinder->fPhiT[i];
+        nchp++;
+     }
+  }
+  printf(" fNtS %i : nchp %i -> %i\n", jetFinder->fNtS, nchp, jetFinder->fNtS - nchp);
+}
+
+void AliEMCALJetMicroDst::FillJetsControl()
+{  // see FillJets(AliEMCALJetFinder* jetFinder) and FillPartons
+  hNJet->Fill(njet);
+  for(Int_t i=0; i<njet; i++){
+     hPtJet->Fill(jet[i]);
+     hEtaPhiJet->Fill(jetaw[i],jphiw[i]); 
+  }
+
+  for(Int_t i=0; i < npart; i++){
+     hEtaPhiPart->Fill(xeta[i], xphi[i]);
+     hPtPart->Fill(xpt[i]);
+  }
+
+  Double_t sum = 0.0;
+  hNcell->Fill(ncell);
+  for(Int_t i=0; i < ncell; i++){
+    hCellId->Fill(idcell[i]);
+    hCellEt->Fill(etcell[i]);
+    sum += Double_t(etcell[i]);
+  }
+  hSumEt->Fill(sum);
+
+  sum = 0.0;
+  hNgrid->Fill(ngrid);
+  for(Int_t i=0; i < ngrid; i++){
+    hGridId->Fill(idgrid[i]);
+    hGridEt->Fill(etgrid[i]);
+    sum += Double_t(etgrid[i]);
+  }
+  hSumEtGrForJF->Fill(sum);
+}
+
 Int_t AliEMCALJetMicroDst::GetEntry(Int_t entry)
 { // Read contents of entry.
    if (!fTree) {
@@ -249,10 +499,10 @@ Int_t AliEMCALJetMicroDst::GetEntry(Int_t entry)
 
 Bool_t AliEMCALJetMicroDst::GetParton(Int_t i, Float_t& pt, Float_t& eta, Float_t& phi)
 {
-  if(i>=0 && i<fNpart) {
-    pt  = fXpt[i];
-    eta = fXeta[i];
-    phi = fXphi[i];
+  if(i>=0 && i<npart) {
+    pt  = xpt[i];
+    eta = xeta[i];
+    phi = xphi[i];
     return kTRUE;
   } else return kFALSE; 
 }
@@ -269,14 +519,14 @@ Bool_t AliEMCALJetMicroDst::GetParton(Int_t i, TVector3& vec)
 
 Bool_t AliEMCALJetMicroDst::GetJet(Int_t i,Int_t mode,Float_t& pt, Float_t& eta, Float_t& phi)
 {// mode=1(W) mode=any(L)
-  if(i>=0 && i<fNjet) {
-    pt  = fJet[i];
+  if(i>=0 && i<njet) {
+    pt  = jet[i];
     if(mode==1) {
-      eta = fJetaW[i];
-      phi = fJphiW[i];
+      eta = jetaw[i];
+      phi = jphiw[i];
     } else {
-      eta = fJetaL[i];
-      phi = fJphiL[i];
+      eta = jetal[i];
+      phi = jphil[i];
     }
     return kTRUE;
   } else return kFALSE; 
@@ -301,20 +551,20 @@ void AliEMCALJetMicroDst::Test()
   for(Int_t i=0; i<nentries; i++){
     nb = fTree->GetEntry(i);  
     nbytes += nb;
-    for(Int_t j=0; j<fNpart; j++){
-      hEtaPhiPart->Fill(fXeta[j], fXphi[j]);
-      hPtPart->Fill(fXpt[j]);
+    for(Int_t j=0; j<npart; j++){
+      hEtaPhiPart->Fill(xeta[j], xphi[j]);
+      hPtPart->Fill(xpt[j]);
     }
 
-    hNJet->Fill(fNjet);
-    if(fNjet){
-      for(Int_t j=0; j<fNjet; j++) {
-        hPtJet->Fill(fJet[j]);
-        hEtaPhiJet->Fill(fJetaW[j],fJphiW[j]); 
+    hNJet->Fill(njet);
+    if(njet){
+      for(Int_t j=0; j<njet; j++) {
+        hPtJet->Fill(jet[j]);
+        hEtaPhiJet->Fill(jetaw[j],jphiw[j]); 
       }
     }
   }
-  printf("\n<I> AliEMCALJetMicroDst::Test() -> Entries %5i Bytes %10i\n", nentries, nb);
+  printf("\n<I> AliEMCALJetMicroDst::Test() -> Entries %5i Bytes %10i\n", nentries, nbytes);
 }
 
 void AliEMCALJetMicroDst::FillVector(Float_t pt, Float_t eta, Float_t phi, TVector3& vec)
@@ -328,6 +578,109 @@ void AliEMCALJetMicroDst::FillVector(Float_t pt, Float_t eta, Float_t phi, TVect
   vec.SetXYZ(px, py, pz);
 }
 
+void AliEMCALJetMicroDst::GetEtaPhi(Int_t id, Double_t &eta, Double_t &phi)
+{ // see AliEMCALGeometry 
+  static Int_t ieta, iphi, nphi=144, neta=96;
+  static Double_t phiMax=TMath::Pi(), phiMin=phiMax/3.;
+  static Double_t phiStep=(phiMax-phiMin)/nphi, phiBeg = phiMin + phiStep/2.; 
+  static Double_t etaMax=0.7, etaMin=-etaMax;
+  static Double_t etaStep=(etaMax-etaMin)/neta, etaBeg = etaMin + etaStep/2.;
+
+  ieta = (id-1)/nphi + 1; // id = nphi*(ieta-1) + iphi
+  iphi = id - nphi*(ieta-1);
+  if(ieta<=0 || ieta>neta) {
+    printf("<E> wrong id %i(ieta %i,iphi %i) : nphi %i neta %i\n", 
+    id,iphi,ieta, nphi,neta);
+    assert(0);
+  }
+
+  eta  = etaBeg + etaStep*(ieta-1);
+  phi  = phiBeg + phiStep*(iphi-1);
+}
+
+TVector3& AliEMCALJetMicroDst::GetCellVector(Int_t i)
+{
+  static Double_t eta,phi;
+  static TVector3 vec;
+  vec.SetXYZ(0.,0.,0.);
+  if(i>=0 && i<ncell) {
+     GetEtaPhi(idcell[i], eta, phi);
+     vec.SetPtEtaPhi(Double_t(etcell[i]),eta,phi);
+  }
+  return vec;
+}
+
+TVector3& AliEMCALJetMicroDst::GetGridVector(Int_t i)
+{
+  static Double_t eta,phi;
+  static TVector3 vec;
+  vec.SetXYZ(0.,0.,0.);
+  if(i>=0 && i<ngrid) {
+     GetEtaPhi(idgrid[i], eta, phi);
+     vec.SetPtEtaPhi(Double_t(etgrid[i]),eta,phi);
+  }
+  return vec;
+}
+
+Double_t AliEMCALJetMicroDst::GetSumInCone(TVector3 &jet,Int_t nc, Float_t *et,Float_t *eta,Float_t *phi, Double_t cellEtCut, Double_t rJet)
+{
+  static Double_t sum=0.;
+  static TVector3 cell(0., 0., 0.);
+  if(nc<=0 || et==0 || eta==0 || phi==0) {
+    cout<<"<E> AliEMCALJetMicroDst::GetSumInCone : nc "
+        <<nc<<" ar "<<et<<eta<<phi<<endl;
+    return -1.;
+  }
+
+  sum=0.;
+  //  jet.SetPtEtaPhi(jet[0],jetaw[0],jphiw[0]); // must be one jet !!
+  printf(" jet.Mag() %f : njet %i\n", jet.Mag(), njet);
+  for(Int_t i=0; i<nc; i++){
+    if(et[i] < cellEtCut)    continue;
+    cell.SetPtEtaPhi(et[i], eta[i], phi[i]);
+    if(jet.DeltaR(cell) > rJet) continue;
+    sum += et[i];
+  }
+  printf(" Sum %f \n", sum);
+  //  assert(0);
+  return sum;
+}
+
+Double_t AliEMCALJetMicroDst::GetEmcalEtInCone(TVector3 &jet, Double_t cellEtCut, Double_t rJet)
+{
+  Int_t nc = ncell;
+  if(nc<=0) return 0.;
+
+  Float_t *et=etcell, *eta=new Float_t[nc], *phi=new Float_t[nc];
+  Double_t etaCell=0., phiCell=0., ET=0;
+
+  for(Int_t i=0; i<nc; i++) {
+    GetEtaPhi(idcell[i], etaCell, phiCell);
+    eta[i] = etaCell;
+    phi[i] = phiCell;
+  }
+
+  ET = GetSumInCone(jet, nc, et,eta,phi, cellEtCut,rJet);
+  delete [] eta;
+  delete [] phi;
+
+  return ET;
+}
+
+Double_t AliEMCALJetMicroDst::GetTpcPtInCone(TVector3 &jet,Double_t cellEtCut, Double_t rJet)
+{
+  if(nchp<=0) return 0.;
+  return GetSumInCone(jet, nchp, ppt,peta,pphi, cellEtCut,rJet);
+}
+
+Double_t AliEMCALJetMicroDst::GetSum(Int_t n, Float_t *ar, Double_t cut)
+{ // 25-apr-2003
+  Double_t sum=0.0;
+  if(n<=0 || ar==0) return sum;
+  for(Int_t i=0; i<n; i++) {if(ar[i]>=cut) sum += Double_t(ar[i]);}
+  return sum;
+}
+
 void AliEMCALJetMicroDst::Close()
 {
   fFile->Write(); 
@@ -344,6 +697,13 @@ void AliEMCALJetMicroDst::Browse(TBrowser* b)
    //   TObject::Browse(b);
 }
 
+Bool_t  AliEMCALJetMicroDst::IsPythiaDst()
+{
+  TString st(GetTitle());
+  if(st.Contains("py",TString::kIgnoreCase)||!st.Contains("hijing", TString::kIgnoreCase)) return kTRUE;
+  else return kFALSE;
+}
+
 Bool_t  AliEMCALJetMicroDst::IsFolder() const
 {
   if(fTree || fListHist) return kTRUE;
index 67d53325cb354693699c61cc5f3f68602e774397..55bf53f97ecec21f078b8d1bbe95f83381ea9a35 100644 (file)
@@ -23,51 +23,86 @@ class TBrowser;
 class AliEMCALJetMicroDst: public TNamed {
 
   private:
-  Int_t  fDebug; 
+  Int_t   fDebug;
 
-  protected:
+  public:
   TFile*  fFile;
   TTree*  fTree;
   TString fName;
   TList*  fListHist;    //!
+  TString fFileName;    // for convenience
 
+  Float_t decone;   //! for EMCAL
+  Float_t ptcone;   //! for ch.particles 
   // For partons after hard scattering
-  Int_t   fNpart;
-  Float_t fXpt[4]; 
-  Float_t fXeta[4]; 
-  Float_t fXphi[4];
+  Int_t   npart;
+  Float_t xpt[4];  //[npart]
+  Float_t xeta[4]; //[npart]
+  Float_t xphi[4]; //[npart]
   // Jet 
-  Int_t   fNjet;
-  Float_t fJet[10]; 
-  Float_t fJetaL[10]; 
-  Float_t fJphiL[10];
-  Float_t fJetaW[10]; 
-  Float_t fJphiW[10];
+  Int_t   njet;
+  Float_t jet[10];   //[njet]
+  Float_t jetal[10]; //[njet]
+  Float_t jphil[10]; //[njet]
+  Float_t jetaw[10]; //[njet]
+  Float_t jphiw[10]; //[njet]
   // Charge particle in jet ??
+  // eT in EMCAL itself - 24-jan-2003
+  Int_t   ncell;         // 96*144 =13824 
+  Int_t   idcell[13824]; //[ncell]
+  Float_t etcell[13824]; //[ncell] : de = det*sf
+  // eT in EMCAL grid for jet finder 
+  Int_t   ngrid;         // 96*144 =13824 
+  Int_t   idgrid[13824]; //[ngrid]
+  Float_t etgrid[13824]; //[ngrid]
+  // charge particle which hit to EMCAL - 28-jan-2003
+  Int_t   nchp;
+  Int_t   pid[20000];  //[nchp]
+  Float_t ppt[20000];  //[nchp]
+  Float_t peta[20000]; //[nchp]
+  Float_t pphi[20000]; //[nchp]
 
   public:
   AliEMCALJetMicroDst(char *name="jetMicroDst",
   char *tit="jet Micro Dst for preparation of proposal");
   virtual ~AliEMCALJetMicroDst();
-  virtual Bool_t  Create(TFile *file);
-  virtual Bool_t  Create(const char  *fname);
-  virtual void    Fill(AliRun *run, AliEMCALJetFinder* jetFinder);
-  virtual void    FillPartons(AliGenHijingEventHeader *header);
-  virtual void    FillPartons();
-  virtual void    FillJets(AliEMCALJetFinder* jetFinder);
+  Bool_t  Create(TFile *file);
+  Bool_t  Create(const char  *fname);
+  void    Fill(AliRun *run=0, AliEMCALJetFinder* jetFinder=0, Int_t modeFilling=0);
+  void    FillPartons(AliGenHijingEventHeader *header);
+  void    FillPartons();
+  void    FillJets(AliEMCALJetFinder* jetFinder);
+  void    FillEtForEMCAL(AliEMCALJetFinder* jetFinder);
+  void    FillEtForGrid(AliEMCALJetFinder* jetFinder);
+  void    FillArrays(TH2* hid, Int_t &n, Int_t *id, Float_t *et);
+  void    FillChargeParticles(AliEMCALJetFinder* jetFinder);  
+  
+  void    FillJetsControl(); // 18-jan-2003
 
-  virtual Bool_t  Open(const char  *fname);    // *MENU* 
-  virtual Bool_t  Initialize(TFile *file);
-  virtual void    Print(Option_t* option="") const;    // *MENU* 
-  virtual Int_t   GetEntry(Int_t entry);
-  virtual void    Test();
-  Int_t   GetNpart() {return fNpart;}
+  Bool_t  Open(const Int_t mode=1) {return Open(DefineName(mode));}  // *MENU* 
+  Bool_t  Open(const char  *fname);                                  // *MENU* 
+  const Char_t* DefineName(const Int_t mode=1);                      // *MENU*
+  Bool_t  Initialize(TFile *file);
+  void    Print(Option_t* option="") const;                          // *MENU* 
+  Int_t   GetEntry(Int_t entry);
+  void    Test();
+  Int_t   GetNpart() {return npart;}
   Bool_t  GetParton(Int_t i, Float_t& pt, Float_t& eta, Float_t& phi);
   Bool_t  GetParton(Int_t i, TVector3& vec);
-  Int_t   GetNjet() {return fNjet;} 
+  Int_t   GetNjet() {return njet;} 
   Bool_t  GetJet(Int_t i,Int_t mode, Float_t& pt,Float_t& eta,Float_t& phi);
   Bool_t  GetJet(Int_t i,Int_t mode, TVector3& vec);
   static  void FillVector(Float_t pt, Float_t eta, Float_t phi, TVector3& vec);
+  void    GetEtaPhi(Int_t id, Double_t &eta, Double_t &phi);
+  TVector3& GetCellVector(Int_t i);
+  TVector3& GetGridVector(Int_t i);
+  // 13-apr-2003
+  Double_t GetSumInCone(TVector3 &jet, Int_t nc, Float_t *et,Float_t *eta,Float_t *phi, Double_t cellEtCut, Double_t rJet);
+  Double_t GetEmcalEtInCone(TVector3 &jet, Double_t cellEtCut=0.0, Double_t rJet=0.5);
+  Double_t GetTpcPtInCone(TVector3 &jet, Double_t cellEtCut=0.0, Double_t rJet=0.5);
+  Double_t GetSum(Int_t n, Float_t *ar, Double_t cut=0.0);
+  Double_t GetSumEmcal(Double_t cut=0.0) {return GetSum(ncell, etcell, cut);}
+  Double_t GetSumTpc(Double_t cut=0.0) {return GetSum(nchp, ppt, cut);}
 
   void    SetDebug(Int_t flag) {fDebug = flag;}
   Float_t GetDebug() const  {return fDebug;}
@@ -76,6 +111,7 @@ class AliEMCALJetMicroDst: public TNamed {
   TFile* GetFile() {return fFile;}
   void   Close();
 
+  Bool_t  IsPythiaDst();
   virtual Bool_t  IsFolder() const;
   virtual void Browse(TBrowser* b);
 
@@ -85,3 +121,7 @@ class AliEMCALJetMicroDst: public TNamed {
 };
 
 #endif // AliEMCALJETMICRODST_H
+/*
+What to do
+1. Common info about event
+ */