+
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-/* $Id$ */
-//*-- Author: Andreas Morsch (CERN)
+/*
+$Log$
+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
+
+*/
+
+//*-- Authors: Andreas Morsch (CERN)
+//* J.L. Klay (LBL)
+//* Aleksei Pavlinov (WSU)
+
+// From root ...
#include <TClonesArray.h>
#include <TTree.h>
+#include <TBranchElement.h>
#include <TFile.h>
#include <TH2.h>
#include <TAxis.h>
#include <TParticle.h>
#include <TParticlePDG.h>
+
+// From AliRoot ...
#include "AliEMCALJetFinder.h"
+#include "AliEMCALFast.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALHit.h"
+#include "AliEMCALDigit.h"
+#include "AliEMCALDigitizer.h"
+#include "AliEMCALHadronCorrection.h"
#include "Ecommon.h"
#include "AliRun.h"
+#include "AliMagF.h"
+#include "AliMagFCM.h"
#include "AliEMCAL.h"
#include "AliHeader.h"
+#include "AliPDG.h"
ClassImp(AliEMCALJetFinder)
AliEMCALJetFinder::AliEMCALJetFinder()
{
// Default constructor
- fJets = 0;
- fNjets = 0;
- fLego = 0;
+ fJets = 0;
+ fNjets = 0;
+ fLego = 0;
+ fTrackList = 0;
+ fPtT = 0;
+ fEtaT = 0;
+ fPhiT = 0;
+ fHadronCorrector = 0;
}
AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
// Constructor
fJets = new TClonesArray("AliEMCALJet",10000);
fNjets = 0;
- for (Int_t i=0; i<30000; i++)
+ for (Int_t i = 0; i < 30000; i++)
{
fEtCell[i] = 0.;
fEtaCell[i] = 0.;
fPhiCell[i] = 0.;
}
fLego = 0;
+ fTrackList = 0;
+ fTrackList = 0;
+ fPtT = 0;
+ fEtaT = 0;
+ fPhiT = 0;
+ fHadronCorrector = 0;
+//
+ SetPtCut();
+ SetMomentumSmearing();
+ SetEfficiencySim();
+ SetDebug();
+ SetHadronCorrection();
+ SetSamplingFraction();
+ SetIncludeK0andN();
}
+
#ifndef WIN32
# define jet_finder_ua1 jet_finder_ua1_
# define hf1 hf1_
jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
min_move, max_move, mode, prec_bg, ierror);
// Write result to output
+ Int_t njet = Njets();
+
+ for (Int_t nj=0; nj<njet; nj++)
+ {
+
+
+ fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
+ JetPhiW(nj),
+ JetEtaW(nj));
+ }
+ FindTracksInJetCone();
WriteJets();
}
// Set grid cell size
EMCALCELLGEO.etaCellSize = eta;
EMCALCELLGEO.phiCellSize = phi;
+ fDeta = eta;
+ fDphi = phi;
}
void AliEMCALJetFinder::SetConeRadius(Float_t par)
{
// Set jet cone radius
EMCALJETPARAM.coneRad = par;
+ fConeRadius = par;
}
void AliEMCALJetFinder::SetEtSeed(Float_t par)
{
// Set et cut for seeds
EMCALJETPARAM.etSeed = par;
+ fEtSeed = par;
}
void AliEMCALJetFinder::SetMinJetEt(Float_t par)
{
// Set minimum jet et
EMCALJETPARAM.ejMin = par;
+ fMinJetEt = par;
}
void AliEMCALJetFinder::SetMinCellEt(Float_t par)
{
// Set et cut per cell
EMCALJETPARAM.etMin = par;
+ fMinCellEt = par;
+}
+
+void AliEMCALJetFinder::SetPtCut(Float_t par)
+{
+// Set pT cut on charged tracks
+ fPtCut = par;
}
const char* file = (pK->GetCurrentFile())->GetName();
// I/O
AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
+
+ if (fDebug > 1)
printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
+
if (fJets && gAlice->TreeR()) {
pEMCAL->MakeBranchInTree(gAlice->TreeR(),
- "Jets",
+ "EMCALJets",
&fJets,
kBufferSize,
file);
}
Int_t njet = Njets();
- for (Int_t nj=0; nj<njet; nj++)
+ for (Int_t nj = 0; nj < njet; nj++)
{
- AliEMCALJet* jet = new AliEMCALJet(JetEnergy(nj),
- JetPhiW(nj),
- JetEtaW(nj));
- AddJet(*jet);
- delete jet;
+ AddJet(*fJetT[nj]);
+ delete fJetT[nj];
}
+
Int_t nev = gAlice->GetHeader()->GetEvent();
gAlice->TreeR()->Fill();
char hname[30];
void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
{
//
-// Fill Cells with hit information
+// Fill Cells with track information
//
//
ResetMap();
//
// Access particle information
Int_t npart = (gAlice->GetHeader())->GetNprimary();
- for (Int_t part=2; part<npart; part++) {
+// Create track list
+//
+// 0: not selected
+// 1: selected for jet finding
+// 2: inside jet
+// ....
+ if (fTrackList) delete[] fTrackList;
+ if (fPtT) delete[] fPtT;
+ if (fEtaT) delete[] fEtaT;
+ if (fPhiT) delete[] fPhiT;
+
+ fTrackList = new Int_t [npart];
+ fPtT = new Float_t[npart];
+ fEtaT = new Float_t[npart];
+ fPhiT = new Float_t[npart];
+ fNt = npart;
+ for (Int_t part = 2; part < npart; part++) {
TParticle *MPart = gAlice->Particle(part);
Int_t mpart = MPart->GetPdgCode();
Int_t child1 = MPart->GetFirstDaughter();
Float_t pT = MPart->Pt();
+ Float_t p = MPart->P();
Float_t phi = MPart->Phi();
- Float_t theta = MPart->Theta();
- Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
-// if (part == 6 || part == 7)
-// {
-// printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f",
-// part-5, pT, eta, phi);
-// }
-
- if (pT == 0.) continue;
+ Float_t eta = MPart->Eta();
+
+ if (fDebug > 0) {
+ if (part == 6 || part == 7)
+ {
+ printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f",
+ part-5, pT, eta, phi);
+ }
+
+// if (mpart == 21)
+
+// printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
+// part, mpart, pT, eta, phi);
+ }
+
+
+ fTrackList[part] = 0;
+ fPtT[part] = pT;
+ fEtaT[part] = eta;
+ fPhiT[part] = phi;
+
+ if (part < 2) continue;
+ if (pT == 0 || pT < fPtCut) continue;
+ TParticlePDG* pdgP = 0;
// charged or neutral
if (ich == 0) {
- TParticlePDG* pdgP = MPart->GetPDG();
- if (pdgP->Charge() == 0) continue;
+ pdgP = MPart->GetPDG();
+ if (pdgP->Charge() == 0) {
+ if (!fK0N) {
+ continue;
+ } else {
+ if (mpart != kNeutron &&
+ mpart != kNeutronBar &&
+ mpart != kK0Long) continue;
+ }
+ }
}
// skip partons
if (TMath::Abs(mpart) <= 6 ||
if (phi*180./TMath::Pi() > 120.) continue;
// final state only
if (child1 >= 0 && child1 < npart) continue;
-// printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
-// part, mpart, child1, eta, phi, pT);
- fLego->Fill(eta, phi, pT);
+//
+ if (fDebug > 1)
+ printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
+ part, mpart, child1, eta, phi, pT);
+//
+//
+// phi propagation
+
+ Bool_t curls = kFALSE;
+ Float_t dphi = PropagatePhi(pT, pdgP->Charge(), curls);
+ if (curls) continue;
+ phi += dphi;
+//
+// Momentum smearing goes here ...
+//
+ if (fSmear) {
+ p = AliEMCALFast::SmearMomentum(1,p);
+ }
+//
+// Tracking Efficiency and TPC acceptance goes here ...
+ if (fEffic) {
+ Float_t eff = AliEMCALFast::Efficiency(1,p);
+ if (AliEMCALFast::RandomReject(eff)) continue;
+ }
+//
+// Correction of Hadronic Energy goes here ...
+//
+ if (fHCorrection) {
+ if (!fHadronCorrector)
+ Fatal("AliEMCALJetFinder",
+ "Hadronic energy correction required but not defined !");
+ Float_t dpH = fHadronCorrector->GetEnergy(p, eta, 7);
+ p -= dpH;
+ }
+//
+// More to do ? Just think about it !
+//
+ fTrackList[part] = 1;
+//
+ fLego->Fill(eta, phi, p);
} // primary loop
DumpLego();
}
void AliEMCALJetFinder::FillFromHits(Int_t flag)
{
//
-// Fill Cells with track information
+// Fill Cells with hit information
//
//
ResetMap();
//
// Access hit information
AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
-// AliEMCALGeometry* geom =
-// AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
TTree *treeH = gAlice->TreeH();
Int_t ntracks = (Int_t) treeH->GetEntries();
Float_t y = mHit->Y(); // y-pos
Float_t z = mHit->Z(); // z-pos
Float_t eloss = mHit->GetEnergy(); // deposited energy
-// Int_t index = mHit->GetId(); // cell index
-// Float_t eta, phi;
-// geom->EtaPhiFromIndex(index, eta, phi);
- Float_t r = TMath::Sqrt(x*x+y*y);
- Float_t theta = TMath::ATan2(r,z);
- Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
- Float_t phi = TMath::ATan2(y,x);
- fLego->Fill(eta, phi, eloss);
-// if (eloss > 1.) printf("\nx,y,z %f %f %f %f %f",
-// r, z, eta, phi, eloss);
-// printf("\n Max %f", fLego->GetMaximum());
+
+ Float_t r = TMath::Sqrt(x*x+y*y);
+ Float_t theta = TMath::ATan2(r,z);
+ Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
+ Float_t phi = TMath::ATan2(y,x);
+ fLego->Fill(eta, phi, fSamplingF*eloss);
} // Hit Loop
} // Track Loop
DumpLego();
}
+void AliEMCALJetFinder::FillFromDigits(Int_t flag)
+{
+//
+// Fill Cells with digit information
+//
+//
+ ResetMap();
+
+ if (!fLego) BookLego();
+ if (flag == 0) fLego->Reset();
+ Int_t nbytes;
+
+
+//
+// Connect digits
+//
+ TClonesArray* digs = new TClonesArray("AliEMCALDigit", 10000);
+ TTree *treeD = gAlice->TreeD();
+ TBranchElement* branchDg = (TBranchElement*)
+ treeD->GetBranch("EMCAL");
+
+ if (!branchDg) Fatal("AliEMCALJetFinder",
+ "Reading digits requested but no digits in file !");
+
+ branchDg->SetAddress(&digs);
+ Int_t nent = (Int_t) branchDg->GetEntries();
+//
+// Connect digitizer
+//
+ AliEMCALDigitizer* digr = new AliEMCALDigitizer();
+ TBranchElement* branchDr = (TBranchElement*)
+ treeD->GetBranch("AliEMCALDigitizer");
+ branchDr->SetAddress(&digr);
+//
+//
+ nbytes += branchDg->GetEntry(0);
+ nbytes += branchDr->GetEntry(0);
+//
+// Get digitizer parameters
+ Float_t towerADCped = digr->GetTowerpedestal();
+ Float_t towerADCcha = digr->GetTowerchannel();
+ Float_t preshoADCped = digr->GetPreShopedestal();
+ Float_t preshoADCcha = digr->GetPreShochannel();
+
+ AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+ AliEMCALGeometry* geom =
+ AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
+
+ if (fDebug) {
+ Int_t ndig = digs->GetEntries();
+ printf("\n Number of Digits: %d %d\n", ndig, nent);
+ printf("\n Parameters: %f %f %f %f\n",
+ towerADCped, towerADCcha, preshoADCped, preshoADCcha );
+ printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
+ }
+
+//
+// Loop over digits
+ AliEMCALDigit* sdg;
+ TIter next(digs);
+ while ((sdg = (AliEMCALDigit*)(next())))
+ {
+ Double_t pedestal;
+ Double_t channel;
+ if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
+ {
+ pedestal = preshoADCped;
+ channel = preshoADCcha;
+ } else {
+ pedestal = towerADCped;
+ channel = towerADCcha;
+ }
+
+ Float_t eta = sdg->GetEta();
+ Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
+ Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
+
+ if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
+ eta, phi, amp, sdg->GetAmp());
+
+ fLego->Fill(eta, phi, fSamplingF*amp);
+ } // digit loop
+//
+// Dump lego hist
+ DumpLego();
+}
+
+
+void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
+{
+//
+// Fill Cells with hit information
+//
+//
+ ResetMap();
+
+ if (!fLego) BookLego();
+// Reset
+ if (flag == 0) fLego->Reset();
+
+//Flag charged tracks passing through TPC which
+//are associated to EMCAL Hits
+ BuildTrackFlagTable();
+
+//
+// Access particle information
+ TTree *treeH = gAlice->TreeH();
+ Int_t ntracks = (Int_t) treeH->GetEntries();
+
+ if (fPtT) delete[] fPtT;
+ if (fEtaT) delete[] fEtaT;
+ if (fPhiT) delete[] fPhiT;
+
+ fPtT = new Float_t[ntracks];
+ fEtaT = new Float_t[ntracks];
+ fPhiT = new Float_t[ntracks];
+ fNt = ntracks;
+
+ for (Int_t track = 0; track < ntracks; track++) {
+ TParticle *MPart = gAlice->Particle(track);
+ Float_t pT = MPart->Pt();
+ Float_t phi = MPart->Phi();
+ Float_t eta = MPart->Eta();
+
+ if(fTrackList[track]) {
+ fPtT[track] = pT;
+ fEtaT[track] = eta;
+ fPhiT[track] = phi;
+
+ if (track < 2) continue; //Colliding particles?
+ if (pT == 0 || pT < fPtCut) continue;
+ fLego->Fill(eta, phi, pT);
+ }
+ } // track loop
+ DumpLego();
+}
+
+void AliEMCALJetFinder::BuildTrackFlagTable() {
+
+// Method to generate a lookup table for TreeK particles
+// which are linked to hits in the EMCAL
+//
+// --Author: J.L. Klay
+//
+ printf("\n Building track flag table.\n");
+
+
+// Access hit information
+ AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
+
+ TTree *TK = gAlice->TreeK(); // Get the number of entries in the kine tree
+ Int_t nKTrks = (Int_t) TK->GetEntries(); // (Number of particles created somewhere)
+
+ if(fTrackList) delete[] fTrackList; //Make sure we get rid of the old one
+ fTrackList = new Int_t[nKTrks]; //before generating a new one
+
+ for (Int_t i = 0; i < nKTrks; i++) { //Initialize members to 0
+ fTrackList[i] = 0;
+ }
+
+ TTree *treeH = gAlice->TreeH();
+ Int_t ntracks = (Int_t) treeH->GetEntries();
+//
+// Loop over tracks
+//
+ Int_t nbytes = 0;
+
+ for (Int_t track=0; track<ntracks;track++) {
+ gAlice->ResetHits();
+ nbytes += treeH->GetEvent(track);
+ if (pEMCAL) {
+
+//
+// Loop over hits
+//
+ for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
+ mHit;
+ mHit=(AliEMCALHit*) pEMCAL->NextHit())
+ {
+ Int_t iTrk = mHit->Track(); // track number
+ 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);
+ TParticlePDG *trkPdg = trkPart->GetPDG();
+ Int_t trkCode = trkPart->GetPdgCode();
+ Double_t trkChg;
+ if (trkCode < 10000) { //Big Ions cause problems for
+ trkChg = trkPdg->Charge(); //this function. Since they aren't
+ } else { //likely to make it very far, set
+ trkChg = 0.0; //their charge to 0 for the Flag test
+ }
+ Float_t vxTrk = trkPart->Vx();
+ Float_t vyTrk = trkPart->Vy();
+ Float_t vrTrk = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
+ fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,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
+ TParticlePDG *ancPdg = ancPart->GetPDG();
+ Int_t ancCode = ancPart->GetPdgCode();
+ Double_t ancChg;
+ if (ancCode < 10000) {
+ ancChg = ancPdg->Charge();
+ } else {
+ ancChg = 0.0;
+ }
+ Float_t vxAnc = ancPart->Vx();
+ Float_t vyAnc = ancPart->Vy();
+ Float_t vrAnc = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
+ fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
+ ancestor = ancPart->GetFirstMother(); //Get the next ancestor
+ }
+
+ //Determine the origin point of the primary particle associated with the hit
+ TParticle *primPart = gAlice->Particle(idprim);
+ TParticlePDG *primPdg = primPart->GetPDG();
+ Int_t primCode = primPart->GetPdgCode();
+ Double_t primChg;
+ if (primCode < 10000) {
+ primChg = primPdg->Charge();
+ } else {
+ primChg = 0.0;
+ }
+ Float_t vxPrim = primPart->Vx();
+ Float_t vyPrim = primPart->Vy();
+ Float_t vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
+ fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
+ } // Hit Loop
+ } //if (pEMCAL)
+ } // Track Loop
+}
+
+Int_t AliEMCALJetFinder
+::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
+
+ Int_t flag = 0;
+ Int_t parton = 0;
+ Int_t neutral = 0;
+
+ if (charge == 0) neutral = 1;
+
+ if (TMath::Abs(code) <= 6 ||
+ code == 21 ||
+ code == 92) parton = 1;
+
+ //It's not a parton, it's charged and it went through the TPC
+ if (!parton && !neutral && radius <= 84.0) flag = 1;
+
+ return flag;
+}
+
+void AliEMCALJetFinder::FindTracksInJetCone()
+{
+//
+// Build list of tracks inside jet cone
+//
+// Loop over jets
+ Int_t njet = Njets();
+ for (Int_t nj = 0; nj < njet; nj++)
+ {
+ Float_t etaj = JetEtaW(nj);
+ Float_t phij = JetPhiW(nj);
+ Int_t nT = 0; // number of associated tracks
+
+// Loop over particles
+ for (Int_t part = 0; part < fNt; part++) {
+ if (!fTrackList[part]) continue;
+ Float_t phi = fPhiT[part];
+ Float_t eta = fEtaT[part];
+ Float_t dr = TMath::Sqrt((etaj-eta)*(etaj-eta) +
+ (phij-phi)*(phij-phi));
+ if (dr < fConeRadius) {
+ fTrackList[part] = nj+2;
+ nT++;
+ } // < ConeRadius ?
+ } // particle loop
+
+ if (nT > 50) nT = 50;
+
+ Float_t* ptT = new Float_t[nT];
+ Float_t* etaT = new Float_t[nT];
+ Float_t* phiT = new Float_t[nT];
+ Int_t iT = 0;
+ Int_t j;
+
+ for (Int_t part = 0; part < fNt; part++) {
+ if (fTrackList[part] == nj+2) {
+ Int_t index = 0;
+ for (j=iT-1; j>=0; j--) {
+ if (fPtT[part] > ptT[j]) {
+ index = j+1;
+ break;
+ }
+ }
+ for (j=iT-1; j>=index; j--) {
+ ptT [j+1] = ptT [j];
+ etaT[j+1] = etaT[j];
+ phiT[j+1] = phiT[j];
+ }
+ ptT [index] = fPtT [part];
+ etaT[index] = fEtaT[part];
+ phiT[index] = fPhiT[part];
+ iT++;
+ } // particle associated
+ } // particle loop
+
+ fJetT[nj]->SetTrackList(nT, ptT, etaT, phiT);
+ delete[] ptT;
+ delete[] etaT;
+ delete[] phiT;
+ } // jet loop loop
+}
+
+Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
+{
+// Propagates phi angle to EMCAL radius
+//
+ Float_t dPhi = 0.;
+// Get field
+ Float_t b = ((AliMagFCM*) gAlice->Field())->SolenoidField();
+// Get EMCAL radius
+ Float_t rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
+//
+//
+// bending radies
+ Float_t rB = 333.56 * pt / b;
+
+//
+// check if particle is curling below EMCAL
+ if (2.*rB < rEMCAL) {
+ curls = kTRUE;
+ return dPhi;
+ }
+//
+// if not calculate delta phi
+ Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
+ dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
+ dPhi = TMath::Sign(dPhi, charge);
+//
+ return dPhi;
+
+}
+
void hf1(Int_t& id, Float_t& x, Float_t& wgt)
{
}
+
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff