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