* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.37 2003/01/23 11:50:04 morsch
-- option for adding energy of all particles (ich == 2)
-- provisions for principle component analysis
-(M. Horner)
-
-Revision 1.36 2003/01/15 19:05:44 morsch
-Updated selection in ReadFromTracks()
-
-Revision 1.35 2003/01/15 04:59:38 morsch
-- TPC eff. from AliEMCALFast
-- Correction in PropagatePhi()
-
-Revision 1.34 2003/01/14 10:50:18 alibrary
-Cleanup of STEER coding conventions
-
-Revision 1.33 2003/01/10 10:48:19 morsch
-SetSamplingFraction() removed from constructor.
-
-Revision 1.32 2003/01/10 10:26:40 morsch
-Sampling fraction initialized from geometry class.
-
-Revision 1.31 2003/01/08 17:13:41 schutz
-added the HCAL section
-
-Revision 1.30 2002/12/09 16:26:28 morsch
-- Nummber of particles per jet increased to 1000
-- Warning removed.
-
-Revision 1.29 2002/11/21 17:01:40 alibrary
-Removing AliMCProcess and AliMC
-
-Revision 1.28 2002/11/20 14:13:16 morsch
-- FindChargedJets() added.
-- Destructor corrected.
-- Geometry cuts taken from AliEMCALGeometry.
-
-Revision 1.27 2002/11/15 17:39:10 morsch
-GetPythiaParticleName removed.
-
-Revision 1.26 2002/10/14 14:55:35 hristov
-Merging the VirtualMC branch to the main development branch (HEAD)
-
-Revision 1.20.4.3 2002/10/10 15:07:49 hristov
-Updating VirtualMC to v3-09-02
-
-Revision 1.25 2002/09/13 10:24:57 morsch
-idem
-
-Revision 1.24 2002/09/13 10:21:13 morsch
-No cast to AliMagFCM.
-
-Revision 1.23 2002/06/27 09:24:26 morsch
-Uncomment the TH1::AddDirectory statement.
-
-Revision 1.22 2002/05/22 13:48:43 morsch
-Pdg code added to track list.
-
-Revision 1.21 2002/04/27 07:43:08 morsch
-Calculation of fDphi corrected (Renan Cabrera)
-
-Revision 1.20 2002/03/12 01:06:23 pavlinov
-Testin output from generator
-
-Revision 1.19 2002/02/27 00:46:33 pavlinov
-Added method FillFromParticles()
-
-Revision 1.18 2002/02/21 08:48:59 morsch
-Correction in FillFromHitFlaggedTrack. (Jennifer Klay)
-
-Revision 1.17 2002/02/14 08:52:53 morsch
-Major updates by Aleksei Pavlinov:
-FillFromPartons, FillFromTracks, jetfinder configuration.
-
-Revision 1.16 2002/02/08 11:43:05 morsch
-SetOutputFileName(..) allows to specify an output file to which the
-reconstructed jets are written. If not set output goes to input file.
-Attention call Init() before processing.
-
-Revision 1.15 2002/02/02 08:37:18 morsch
-Formula for rB corrected.
-
-Revision 1.14 2002/02/01 08:55:30 morsch
-Fill map with Et and pT.
-
-Revision 1.13 2002/01/31 09:37:36 morsch
-Geometry parameters in constructor and call of SetCellSize()
-
-Revision 1.12 2002/01/23 13:40:23 morsch
-Fastidious debug print statement removed.
-
-Revision 1.11 2002/01/22 17:25:47 morsch
-Some corrections for event mixing and bg event handling.
-
-Revision 1.10 2002/01/22 10:31:44 morsch
-Some correction for bg mixing.
-
-Revision 1.9 2002/01/21 16:28:42 morsch
-Correct sign of dphi.
-
-Revision 1.8 2002/01/21 16:05:31 morsch
-- different phi-bin for hadron correction
-- provisions for background mixing.
-
-Revision 1.7 2002/01/21 15:47:18 morsch
-Bending radius correctly in cm.
-
-Revision 1.6 2002/01/21 12:53:50 morsch
-authors
-
-Revision 1.5 2002/01/21 12:47:47 morsch
-Possibility to include K0long and neutrons.
-
-Revision 1.4 2002/01/21 11:03:21 morsch
-Phi propagation introduced in FillFromTracks.
-
-Revision 1.3 2002/01/18 05:07:56 morsch
-- hadronic correction
-- filling of digits
-- track selection upon EMCAL information
-
-*/
+/* $Id$ */
//*-- Authors: Andreas Morsch (CERN)
//* J.L. Klay (LBL)
#include <TROOT.h>
#include <TStyle.h>
#include <TTree.h>
+#include <TBrowser.h>
// From AliRoot ...
#include "AliEMCAL.h"
#include "AliMagFCM.h"
#include "AliRun.h"
#include "AliGenerator.h"
+#include "AliEMCALGetter.h"
// Interface to FORTRAN
#include "Ecommon.h"
+#include "AliMC.h"
ClassImp(AliEMCALJetFinder)
//
//
// Geometry
- AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
- AliEMCALGeometry* geom =
- AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+ //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+ // AliEMCALGeometry* geom =
+ // AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
+ AliEMCALGeometry* geom = gime->EMCALGeometry() ;
// SetSamplingFraction(geom->GetSampling());
}
// I/O
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
if (!fOutFileName) {
//
// output written to input file
//
- AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
- TTree* pK = gAlice->TreeK();
- file = (pK->GetCurrentFile())->GetName();
+ AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
+ TTree* pK = gAlice->TreeK();
+ file = (pK->GetCurrentFile())->GetName();
+ TBranch * jetBranch ;
if (fDebug > 1)
printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
- if (fJets && gAlice->TreeR()) {
- pEMCAL->MakeBranchInTree(gAlice->TreeR(),
- "EMCALJets",
- &fJets,
- kBufferSize,
- file);
+ //if (fJets && gAlice->TreeR()) {
+ if (fJets && gime->TreeR()) {
+ // pEMCAL->MakeBranchInTree(gAlice->TreeR(),
+ jetBranch = gime->TreeR()->Branch("EMCALJets", &fJets, kBufferSize, 0) ;
+ //pEMCAL->MakeBranchInTree(gime->TreeR(),
+ // "EMCALJets",
+ // &fJets,
+ // kBufferSize,
+ // file);
+
+ //Int_t nev = gAlice->GetHeader()->GetEvent();
+ //gAlice->TreeR()->Fill();
+ jetBranch->Fill();
+ //char hname[30];
+ //sprintf(hname,"TreeR%d", nev);
+ //gAlice->TreeR()->Write(hname);
+ //gAlice->TreeR()->Reset();
+ gime->WriteRecPoints("OVERWRITE");
}
- Int_t nev = gAlice->GetHeader()->GetEvent();
- gAlice->TreeR()->Fill();
- char hname[30];
- sprintf(hname,"TreeR%d", nev);
- gAlice->TreeR()->Write(hname);
- gAlice->TreeR()->Reset();
} else {
//
// Output written to user specified output file
//
- TTree* pK = gAlice->TreeK();
- fInFile = pK->GetCurrentFile();
-
- fOutFile->cd();
- char hname[30];
- sprintf(hname,"TreeR%d", fEvent);
- TTree* treeJ = new TTree(hname, "EMCALJets");
- treeJ->Branch("EMCALJets", &fJets, kBufferSize);
- treeJ->Fill();
- treeJ->Write(hname);
- fInFile->cd();
+ //TTree* pK = gAlice->TreeK();
+ TTree* pK = gAlice->TreeK();
+ fInFile = pK->GetCurrentFile();
+
+ fOutFile->cd();
+ char hname[30];
+ sprintf(hname,"TreeR%d", fEvent);
+ TTree* treeJ = new TTree(hname, "EMCALJets");
+ treeJ->Branch("EMCALJets", &fJets, kBufferSize);
+ treeJ->Fill();
+ treeJ->Write(hname);
+ fInFile->cd();
}
ResetJets();
}
void AliEMCALJetFinder::BookLego()
{
//
-// Book histo for discretisation
+// Book histo for discretization
//
//
// Don't add histos to the current directory
if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
- TH2::AddDirectory(0);
- TH1::AddDirectory(0);
+ // TH2::AddDirectory(0); // hists wil be put to the list from gROOT
+ // TH1::AddDirectory(0);
gROOT->cd();
//
// Signal map
fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
"Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
- //! first canvas for drawing
- fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kTRUE);
+ fhSinTheta = new TH1F("fhSinTheta","sin(theta)", fNbinEta, fEtaMin, fEtaMax);
+ TAxis *xax = fhSinTheta->GetXaxis();
+ for(Int_t i=1; i<=fNbinEta; i++) {
+ Double_t eta = xax->GetBinCenter(i);
+ fhSinTheta->Fill(eta, 1./TMath::CosH(eta)); // cosh(eta) = 1./sin(theta)
+ }
+
+ //! first canvas for drawing
+ fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kFALSE);
}
void AliEMCALJetFinder::DumpLego()
// Access particle information
Int_t npart = (gAlice->GetHeader())->GetNprimary();
Int_t ntr = (gAlice->GetHeader())->GetNtrack();
- printf(" : #primary particles %i # tracks %i \n", npart, ntr);
+ printf(" : #primary particles %i # tracks %i : (before) Sum.Et %f\n",
+ npart, ntr, fLego->Integral());
// Create track list
//
// this is for Pythia ??
for (Int_t part = 0; part < npart; part++) {
- TParticle *MPart = gAlice->Particle(part);
+ TParticle *MPart = gAlice->GetMCApp()->Particle(part);
Int_t mpart = MPart->GetPdgCode();
Int_t child1 = MPart->GetFirstDaughter();
Float_t pT = MPart->Pt();
Float_t eta = -100.;
if(pT > 0.001) eta = MPart->Eta();
Float_t theta = MPart->Theta();
- if (fDebug>=2) {
+ if (fDebug>=2 && MPart->GetStatusCode()==1) {
printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
}
if (etc > fMinCellEt) etsum += etc;
}
- printf("\nFillFromTracks: Sum above threshold %f %f \n \n", fMinCellEt, etsum);
+ printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral());
+ printf(" Track selected(fNtS) %i \n", fNtS);
DumpLego();
}
//
// Access hit information
AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
- TTree *treeH = gAlice->TreeH();
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
+ TTree *treeH = gime->TreeH();
Int_t ntracks = (Int_t) treeH->GetEntries();
//
// Loop over tracks
//
Int_t nbytes = 0;
- Double_t etH = 0.0;
+ // Double_t etH = 0.0;
for (Int_t track=0; track<ntracks;track++) {
gAlice->ResetHits();
Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
Float_t phi = TMath::ATan2(y,x);
- if (fDebug >= 11) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
-// printf("\n Hit %f %f %f %f", x, y, z, eloss);
+ if (fDebug >= 21) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
- etH = fSamplingF*eloss*TMath::Sin(theta);
- fLego->Fill(eta, phi, etH);
- // fhLegoEMCAL->Fill(eta, phi, etH);
+ // etH = fSamplingF*eloss*TMath::Sin(theta);
+ fLego->Fill(eta, phi, eloss);
} // Hit Loop
} // Track Loop
- // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
- Float_t etsum = 0;
-
- for(Int_t i=0; i<fLego->GetSize(); i++) {
- (*fhLegoEMCAL)[i] = (*fLego)[i];
- Float_t etc = (*fLego)[i];
- if (etc > fMinCellEt) etsum += etc;
+
+ // Transition from deposit energy to eT (eT = de*SF*sin(theta))
+ Double_t etsum = 0;
+ for(Int_t i=1; i<=fLego->GetNbinsX(); i++){ // eta
+ Double_t sinTheta = fhSinTheta->GetBinContent(i), eT=0;
+ for(Int_t j=1; j<=fLego->GetNbinsY(); j++){ // phi
+ eT = fLego->GetBinContent(i,j)*fSamplingF*sinTheta;
+ fLego->SetBinContent(i,j,eT);
+ // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
+ fhLegoEMCAL->SetBinContent(i,j,eT);
+ if (eT > fMinCellEt) etsum += eT;
+ }
}
+
+ // for(Int_t i=0; i<fLego->GetSize(); i++) {
+ // (*fhLegoEMCAL)[i] = (*fLego)[i];
+ // Float_t etc = (*fLego)[i];
+ // if (etc > fMinCellEt) etsum += etc;
+ // }
- printf("\nFillFromHits: Sum above threshold %f %f \n \n", fMinCellEt, etsum);
-
+ printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum);
// DumpLego(); ??
-
}
void AliEMCALJetFinder::FillFromDigits(Int_t flag)
if (!fLego) BookLego();
if (flag == 0) fLego->Reset();
- Int_t nbytes;
+ Int_t nbytes=0;
//
// Get digitizer parameters
Float_t preADCped = digr->GetPREpedestal();
Float_t preADCcha = digr->GetPREchannel();
- Float_t ecADCped = digr->GetECpedestal();
- Float_t ecADCcha = digr->GetECchannel();
- Float_t hcADCped = digr->GetHCpedestal();
- Float_t hcADCcha = digr->GetHCchannel();
+ Float_t ecADCped = digr->GetECApedestal();
+ Float_t ecADCcha = digr->GetECAchannel();
+ Float_t hcADCped = digr->GetHCApedestal();
+ Float_t hcADCcha = digr->GetHCAchannel();
AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
AliEMCALGeometry* geom =
pedestal = preADCped;
channel = preADCcha;
}
- else if (geom->IsInECAL(sdg->GetId())) {
+ else if (geom->IsInECA(sdg->GetId())) {
pedestal = ecADCped;
channel = ecADCcha;
}
- else if (geom->IsInHCAL(sdg->GetId())) {
+ else if (geom->IsInHCA(sdg->GetId())) {
pedestal = hcADCped;
channel = hcADCcha;
}
fNtS = 0;
for (Int_t track = 0; track < ntracks; track++) {
- TParticle *MPart = gAlice->Particle(track);
+ TParticle *MPart = gAlice->GetMCApp()->Particle(track);
Float_t pT = MPart->Pt();
Float_t phi = MPart->Phi();
Float_t eta = MPart->Eta();
fTrackList[part] = 0;
- MPart = gAlice->Particle(part);
+ MPart = gAlice->GetMCApp()->Particle(part);
mpart = MPart->GetPdgCode();
child1 = MPart->GetFirstDaughter();
child2 = MPart->GetLastDaughter();
// Go through the partons
Int_t statusCode=0;
for (Int_t part = 8; part < npart; part++) {
- TParticle *MPart = gAlice->Particle(part);
+ TParticle *MPart = gAlice->GetMCApp()->Particle(part);
Int_t mpart = MPart->GetPdgCode();
// Int_t child1 = MPart->GetFirstDaughter();
Float_t pT = MPart->Pt();
fTrackList[i] = 0;
}
- TTree *treeH = gAlice->TreeH();
- Int_t ntracks = (Int_t) treeH->GetEntries();
+ AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
+ // TTree *treeH = gAlice->TreeH();
+ TTree *treeH = gime->TreeH();
+ Int_t ntracks = (Int_t) treeH->GetEntries();
//
// Loop over tracks
//
Int_t idprim = mHit->GetPrimary(); // primary particle
//Determine the origin point of this particle - it made a hit in the EMCAL
- TParticle *trkPart = gAlice->Particle(iTrk);
+ TParticle *trkPart = gAlice->GetMCApp()->Particle(iTrk);
TParticlePDG *trkPdg = trkPart->GetPDG();
Int_t trkCode = trkPart->GetPdgCode();
Double_t trkChg;
//Loop through the ancestry of the EMCAL entrance particles
Int_t ancestor = trkPart->GetFirstMother(); //Get track's Mother
while (ancestor != -1) {
- TParticle *ancPart = gAlice->Particle(ancestor); //get particle info on ancestor
+ TParticle *ancPart = gAlice->GetMCApp()->Particle(ancestor); //get particle info on ancestor
TParticlePDG *ancPdg = ancPart->GetPDG();
Int_t ancCode = ancPart->GetPdgCode();
Double_t ancChg;
}
//Determine the origin point of the primary particle associated with the hit
- TParticle *primPart = gAlice->Particle(idprim);
+ TParticle *primPart = gAlice->GetMCApp()->Particle(idprim);
TParticlePDG *primPdg = primPart->GetPDG();
Int_t primCode = primPart->GetPdgCode();
Double_t primChg;
-void AliEMCALJetFinder::SaveBackgroundEvent()
+void AliEMCALJetFinder::SaveBackgroundEvent(Char_t *name)
{
-// Saves the eta-phi lego and the tracklist
+// Saves the eta-phi lego and the tracklist and name of file with BG events
//
if (fLegoB) {
fLegoB->Reset();
(*fLegoB) = (*fLegoB) + (*fLego);
- if(fDebug)
+ // if(fDebug)
printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
fLegoB->Integral(), fLego->Integral());
}
fNtB++;
}
fBackground = 1;
- printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
+ printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
+
+ if(strlen(name) == 0) {
+ TSeqCollection *li = gROOT->GetListOfFiles();
+ TString nf;
+ for(Int_t i=0; i<li->GetSize(); i++) {
+ nf = ((TFile*)li->At(i))->GetName();
+ if(nf.Contains("backgorund")) break;
+ }
+ fBGFileName = nf;
+ } else {
+ fBGFileName = name;
+ }
+ printf("BG file name is \n %s\n", fBGFileName.Data());
}
void AliEMCALJetFinder::InitFromBackground()
fLego->Reset();
(*fLego) = (*fLego) + (*fLegoB);
if(fDebug)
- printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
+ printf("\n AliEMCALJetFinder::InitBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
fLego->Integral(), fLegoB->Integral());
} else {
printf(" => fLego undefined \n");
}
void AliEMCALJetFinder::DrawLego(Char_t *opt)
-{fLego->Draw(opt);}
+{if(fLego) fLego->Draw(opt);}
+
+void AliEMCALJetFinder::DrawLegoBackground(Char_t *opt)
+{if(fLegoB) fLegoB->Draw(opt);}
void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt)
-{fhLegoEMCAL->Draw(opt);}
+{if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt);}
void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
{
if(file==0) file = stdout;
}
fprintf(file,"==== Filling lego ==== \n");
+ fprintf(file,"Sampling fraction %6.3f ", fSamplingF);
fprintf(file,"Smearing %6i ", fSmear);
fprintf(file,"Efficiency %6i\n", fEffic);
fprintf(file,"Hadr.Correct. %6i ", fHCorrection);
fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
} else
fprintf(file,"==== No Bg subtraction ==== \n");
+ //
+ printf("BG file name is %s \n", fBGFileName.Data());
if(file != stdout) fclose(file);
}
if (fWrite) WriteJets();
fEvent++;
}
+// 16-jan-2003 - just for convenience
+void AliEMCALJetFinder::Browse(TBrowser* b)
+{
+ if(fHistsList) b->Add((TObject*)fHistsList);
+}
+
+Bool_t AliEMCALJetFinder::IsFolder() const
+{
+ if(fHistsList) return kTRUE;
+ else return kFALSE;
+}
+
+const Char_t* AliEMCALJetFinder::GetNameOfVariant()
+{// generate the literal string with info about jet finder
+ Char_t name[200];
+ sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f",
+ fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF);
+ TString nt(name);
+ nt.ReplaceAll(" ","");
+ if(fBGFileName.Length()) {
+ Int_t i1 = fBGFileName.Index("kBackground");
+ Int_t i2 = fBGFileName.Index("/0000") - 1;
+ if(i1>=0 && i2>=0) {
+ TString bg(fBGFileName(i1,i2-i1+1));
+ nt += bg;
+ }
+ }
+ printf("<I> Name of variant %s \n", nt.Data());
+ return nt.Data();
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff