X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALJetFinder.cxx;h=b5c6718f5a754cf828281091f27e38b12daee7af;hb=70b76f8245fcad223e5a5ce5cea9e620a4e04f15;hp=0efbcd518fac9b1ab14cc3a77ae3527b5df234e0;hpb=5649a7781bbbdc9c4f72a14a9ac990341d8969b2;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALJetFinder.cxx b/EMCAL/AliEMCALJetFinder.cxx index 0efbcd518fa..b5c6718f5a7 100644 --- a/EMCAL/AliEMCALJetFinder.cxx +++ b/EMCAL/AliEMCALJetFinder.cxx @@ -13,151 +13,56 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -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) //* Aleksei Pavlinov (WSU) +// +// +// #include // From root ... -#include -#include -#include +#include +#include #include +#include +#include #include #include #include -#include -#include #include +#include #include -#include -#include -#include #include #include +#include #include +#include +#include +#include +#include // From AliRoot ... -#include "AliEMCALJetFinder.h" -#include "AliEMCALFast.h" -#include "AliEMCALGeometry.h" -#include "AliEMCALHit.h" +#include "AliEMCAL.h" #include "AliEMCALDigit.h" #include "AliEMCALDigitizer.h" +#include "AliEMCALFast.h" +#include "AliEMCALGeometry.h" #include "AliEMCALHadronCorrection.h" +#include "AliEMCALHit.h" +#include "AliEMCALJetFinder.h" #include "AliEMCALJetMicroDst.h" -#include "AliRun.h" +#include "AliHeader.h" #include "AliMagF.h" #include "AliMagFCM.h" -#include "AliEMCAL.h" -#include "AliHeader.h" -#include "AliPDG.h" - +#include "AliRun.h" +#include "AliGenerator.h" +#include "AliEMCALGetter.h" // Interface to FORTRAN #include "Ecommon.h" +#include "AliMC.h" ClassImp(AliEMCALJetFinder) @@ -192,6 +97,8 @@ AliEMCALJetFinder::AliEMCALJetFinder() fInFile = 0; fEvent = 0; + fRandomBg = 0; + SetParametersForBgSubtraction(); } @@ -238,9 +145,9 @@ AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title) SetEfficiencySim(); SetDebug(); SetHadronCorrection(); - SetSamplingFraction(); SetIncludeK0andN(); + fRandomBg = 0; SetParametersForBgSubtraction(); } @@ -322,11 +229,13 @@ void AliEMCALJetFinder::Init() // // // 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()); +// SetSamplingFraction(geom->GetSampling()); fNbinEta = geom->GetNZ(); fNbinPhi = geom->GetNPhi(); @@ -337,11 +246,15 @@ void AliEMCALJetFinder::Init() fDphi = (fPhiMax-fPhiMin)/fNbinPhi; fDeta = (fEtaMax-fEtaMin)/fNbinEta; fNtot = fNbinPhi*fNbinEta; + fWeightingMethod = kFALSE; + // SetCellSize(fDeta, fDphi); // // I/O if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate"); +// +// } void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000], @@ -368,7 +281,6 @@ void AliEMCALJetFinder::Find() Int_t mode = fMode; Float_t prec_bg = fPrecBg; Int_t ierror; - ResetJets(); // 4-feb-2002 by PAI jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell, @@ -379,10 +291,23 @@ void AliEMCALJetFinder::Find() for (Int_t nj=0; njSetIsWeightedEnergy(fWeightingMethod); + fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) ); + fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj) )); + fJetT[nj]->SetHCEnergy(HCConeEnergy( JetEtaW(nj),JetPhiW(nj) )); + } FindTracksInJetCone(); @@ -390,7 +315,106 @@ void AliEMCALJetFinder::Find() fEvent++; } -Int_t AliEMCALJetFinder::Njets() + +Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi) +{ +// Calculate the energy in the cone +Float_t newenergy = 0.0; +Float_t bineta,binphi; +TAxis *x = fhLegoEMCAL->GetXaxis(); +TAxis *y = fhLegoEMCAL->GetYaxis(); +for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord + { + for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord + { + bineta = x->GetBinCenter(i); + binphi = y->GetBinCenter(j); + if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius) + { + newenergy += fhLegoEMCAL->GetBinContent(i,j); + } + } +} +return newenergy; +} + +Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi) +{ +// Calculate the track energy in the cone +Float_t newenergy = 0.0; +Float_t bineta,binphi; +TAxis *x = fhLegoTracks->GetXaxis(); +TAxis *y = fhLegoTracks->GetYaxis(); +for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord + { + for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord + { + bineta = x->GetBinCenter(i); + binphi = y->GetBinCenter(j); + if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius) + { + newenergy += fhLegoTracks->GetBinContent(i,j); + } + } +} +return newenergy; +} + +Float_t AliEMCALJetFinder::HCConeEnergy(Float_t eta, Float_t phi) +{ +//Float_t newenergy = 0.0; +//Float_t bineta,binphi; +//TAxis *x = fhLegoTracks->GetXaxis(); +//TAxis *y = fhLegoTracks->GetYaxis(); +//for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord +// { +// for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord +// { +// bineta = x->GetBinCenter(i); +// binphi = y->GetBinCenter(j); +// if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius) +// { +// newenergy += fhLegoTracks->GetBinContent(i,j); +// } +// } +//} +//return newenergy; + +return eta*phi*0.0; + +} + + + +Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi) +{ +// Calculate the weighted jet energy + +Float_t newenergy = 0.0; +Float_t bineta,binphi; +TAxis *x = fhLegoEMCAL->GetXaxis(); +TAxis *y = fhLegoEMCAL->GetYaxis(); + + +for (Int_t i = 0 ; i < fNbinEta ; i++) // x coord +{ + for (Int_t j = 0 ; j < fNbinPhi ; j++) // y coord + { + bineta = x->GetBinCenter(i); + binphi = y->GetBinCenter(j); + if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius) + { + newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j); + } + } +} + +return newenergy; + +} + + +Int_t AliEMCALJetFinder::Njets() const { // Get number of reconstructed jets return EMCALJETS.njet; @@ -532,43 +556,52 @@ void AliEMCALJetFinder::WriteJets() } // 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(); } @@ -576,15 +609,15 @@ void AliEMCALJetFinder::WriteJets() 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 @@ -626,8 +659,15 @@ void AliEMCALJetFinder::BookLego() 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() @@ -636,29 +676,37 @@ void AliEMCALJetFinder::DumpLego() // Dump lego histo into array // fNcell = 0; - TAxis* Xaxis = fLego->GetXaxis(); - TAxis* Yaxis = fLego->GetYaxis(); + TAxis* xaxis = fLego->GetXaxis(); + TAxis* yaxis = fLego->GetYaxis(); // fhCellEt->Clear(); Float_t e, eH; for (Int_t i = 1; i <= fNbinEta; i++) { for (Int_t j = 1; j <= fNbinPhi; j++) { + e = fLego->GetBinContent(i,j); - if (e > 0.0) { - Float_t eta = Xaxis->GetBinCenter(i); - Float_t phi = Yaxis->GetBinCenter(j); - fEtCell[fNcell] = e; - fEtaCell[fNcell] = eta; - fPhiCell[fNcell] = phi; - fNcell++; - fhCellEt->Fill(e); - } + if (fRandomBg) { + if (gRandom->Rndm() < 0.5) { + Float_t ebg = 0.28 * TMath::Abs(gRandom->Gaus(0.,1.)); + e += ebg; + } + } + if (e > 0.0) e -= fMinCellEt; + if (e < 0.0) e = 0.; + Float_t eta = xaxis->GetBinCenter(i); + Float_t phi = yaxis->GetBinCenter(j); + fEtCell[fNcell] = e; + fEtaCell[fNcell] = eta; + fPhiCell[fNcell] = phi; + fNcell++; + fhCellEt->Fill(e); if(fhLegoEMCAL) { eH = fhLegoEMCAL->GetBinContent(i,j); if(eH > 0.0) fhCellEMCALEt->Fill(eH); } } // phi } // eta - fNcell--; +// today +// fNcell--; } void AliEMCALJetFinder::ResetMap() @@ -677,6 +725,23 @@ void AliEMCALJetFinder::ResetMap() void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) { +// Which generator +// + const char* name = gAlice->Generator()->GetName(); + enum {kPythia, kHijing, kHijingPara}; + Int_t genType = 0; + + if (!strcmp(name, "Hijing")){ + genType = kHijing; + } else if (!strcmp(name, "Pythia")) { + genType = kPythia; + } else if (!strcmp(name, "HIJINGpara") ||!strcmp(name, "HIGINGpara")) { + genType = kHijingPara; + } + if (fDebug>=2) + printf("Fill tracks generated by %s %d\n", name, genType); + + // // Fill Cells with track information // @@ -693,7 +758,8 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) // 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 // @@ -720,7 +786,7 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) // 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(); @@ -729,22 +795,17 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) 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 (fDebug >= 2) { + if (fDebug >= 2 && genType == kPythia) { if (part == 6 || part == 7) { printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n", 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; @@ -753,12 +814,15 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) fPhiT[part] = phi; fPdgT[part] = mpart; +// final state particles only - if (part < 2) continue; - - // move to fLego->Fill because hadron correction must apply - // if particle hit to EMCAL - 11-feb-2002 - // if (pT == 0 || pT < fPtCut) continue; + if (genType == kPythia) { + if (MPart->GetStatusCode() != 1) continue; + } else if (genType == kHijing) { + if (child1 >= 0 && child1 < npart) continue; + } + + TParticlePDG* pdgP = 0; // charged or neutral pdgP = MPart->GetPDG(); @@ -774,13 +838,12 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) mpart != kK0Long) continue; } } + } else if (ich == 2) { + if (mpart == kNeutron || + mpart == kNeutronBar || + mpart == kK0Long) continue; } -// skip partons - if (TMath::Abs(mpart) <= 6 || - mpart == 21 || - mpart == 92) continue; - if (TMath::Abs(eta)<=0.9) fNChTpc++; // final state only if (child1 >= 0 && child1 < npart) continue; @@ -799,8 +862,8 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) Float_t pw; if (fSmear && TMath::Abs(chTmp)) { pw = AliEMCALFast::SmearMomentum(1,p); - // p changed - take into account when calculate pT, - // pz and so on .. ? + // p changed - take into account when calculate pT, + // pz and so on .. ? pT = (pw/p) * pT; if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw); p = pw; @@ -809,12 +872,11 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) // Tracking Efficiency and TPC acceptance goes here ... Float_t eff; if (fEffic && TMath::Abs(chTmp)) { - // eff = AliEMCALFast::Efficiency(1,p); - eff = 0.95; // for testing 25-feb-2002 + eff = 0.9; // AliEMCALFast::Efficiency(2,p); if(fhEff) fhEff->Fill(p, eff); if (AliEMCALFast::RandomReject(eff)) { - if(fDebug >= 5) printf(" reject due to unefficiency "); - continue; + if(fDebug >= 5) printf(" reject due to unefficiency "); + continue; } } // @@ -826,28 +888,33 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) Bool_t curls = kFALSE; // hit two the EMCAL (no curl) Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0; if(TMath::Abs(chTmp)) { - // hadr. correction only for charge particle - dphi = PropagatePhi(pT, chTmp, curls); - phiHC = phi + dphi; - if (fDebug >= 6) { - printf("\n Delta phi %f pT %f ", dphi, pT); - if (curls) printf("\n !! Track is curling"); - } - if(!curls) fhTrackPtBcut->Fill(pT); - - if (fHCorrection && !curls) { - if (!fHadronCorrector) - Fatal("AliEMCALJetFinder", - "Hadronic energy correction required but not defined !"); - - dpH = fHadronCorrector->GetEnergy(p, eta, 7); - eTdpH = dpH*TMath::Sin(theta); - - if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n", - phi, phiHC, -eTdpH); // correction is negative - fLego->Fill(eta, phiHC, -eTdpH); - fhLegoHadrCorr->Fill(eta, phiHC, eTdpH); - } + // hadr. correction only for charge particle + dphi = PropagatePhi(pT, chTmp, curls); + phiHC = phi + dphi; + if (fDebug >= 6) { + printf("\n Delta phi %f pT %f ", dphi, pT); + if (curls) printf("\n !! Track is curling"); + } + if(!curls) fhTrackPtBcut->Fill(pT); + + if (fHCorrection && !curls) { + if (!fHadronCorrector) + Fatal("AliEMCALJetFinder", + "Hadronic energy correction required but not defined !"); + + dpH = fHadronCorrector->GetEnergy(p, eta, 7); + eTdpH = dpH*TMath::Sin(theta); + + if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n", + phi, phiHC, -eTdpH); // correction is negative + Int_t xbin,ybin; + xbin = fLego->GetXaxis()->FindBin(eta); + ybin = fLego->GetYaxis()->FindBin(phiHC); + cout <<"Hadron Correction affected bin - contents before correction : "<GetBinContent(xbin,ybin)<Fill(eta, phi, -fSamplingF*eTdpH ); + cout <<"Hadron Correction affected bin - contents after correction : "<GetBinContent(xbin,ybin)<Fill(eta, phi, fSamplingF*eTdpH); + } } // // More to do ? Just think about it ! @@ -855,24 +922,32 @@ void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich) if (phi > fPhiMax || phi < fPhiMin ) continue; if(TMath::Abs(chTmp) ) { // charge particle - if (pT > fPtCut && !curls) { - if (fDebug >= 8) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n", - eta , phi, pT, fNtS); - fLego->Fill(eta, phi, pT); - fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking - fTrackList[part] = 1; - fNtS++; - } - } else if(ich==0 && fK0N) { - // case of n, nbar and K0L - if (fDebug >= 9) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n", - eta , phi, pT, fNtS); + if (pT > fPtCut && !curls) { + if (fDebug >= 8) printf("Charge : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n", + eta , phi, pT, fNtS); + fLego->Fill(eta, phi, pT); + fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking + fTrackList[part] = 1; + fNtS++; + } + } else if(ich > 0 || fK0N) { + // case of n, nbar and K0L + if (fDebug >= 9) printf("Neutral : fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n", + eta , phi, pT, fNtS); fLego->Fill(eta, phi, pT); fTrackList[part] = 1; fNtS++; } - } // primary loop + Float_t etsum = 0.; + for(Int_t i=0; iGetSize(); i++) { + Float_t etc = (*fLego)[i]; + if (etc > fMinCellEt) etsum += etc; + } + + printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral()); + printf(" Track selected(fNtS) %i \n", fNtS); + DumpLego(); } @@ -899,13 +974,14 @@ void AliEMCALJetFinder::FillFromHits(Int_t flag) // // 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; trackResetHits(); @@ -927,16 +1003,33 @@ void AliEMCALJetFinder::FillFromHits(Int_t flag) 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", x, y, z, eloss); -// 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) - for(Int_t i=0; iGetSize(); i++) (*fhLegoEMCAL)[i] = (*fLego)[i]; + + // 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; iGetSize(); i++) { + // (*fhLegoEMCAL)[i] = (*fLego)[i]; + // Float_t etc = (*fLego)[i]; + // if (etc > fMinCellEt) etsum += etc; + // } + + printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum); // DumpLego(); ?? } @@ -950,7 +1043,7 @@ void AliEMCALJetFinder::FillFromDigits(Int_t flag) if (!fLego) BookLego(); if (flag == 0) fLego->Reset(); - Int_t nbytes; + Int_t nbytes=0; // @@ -981,10 +1074,10 @@ void AliEMCALJetFinder::FillFromDigits(Int_t flag) // 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 = @@ -1008,11 +1101,11 @@ void AliEMCALJetFinder::FillFromDigits(Int_t flag) 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; } @@ -1070,16 +1163,16 @@ void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag) fNtS = 0; 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(); + TParticle *mPart = gAlice->GetMCApp()->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; - fPdgT[track] = MPart->GetPdgCode(); + fPdgT[track] = mPart->GetPdgCode(); if (track < 2) continue; //Colliding particles? if (pT == 0 || pT < fPtCut) continue; @@ -1095,7 +1188,7 @@ void AliEMCALJetFinder::FillFromParticles() // 26-feb-2002 PAI - for checking all chain // Work on particles level; accept all particle (not neutrino ) - Double_t PX=0, PY=0, PZ=0, E=0; // checking conservation law + Double_t pX=0, pY=0, pZ=0, energy=0; // checking conservation law fNChTpc = 0; ResetMap(); @@ -1115,23 +1208,23 @@ void AliEMCALJetFinder::FillFromParticles() // Go through the particles Int_t mpart, child1, child2, geantPdg; Float_t pT, phi, eta, e=0, px=0, py=0, pz=0; - TParticle *MPart=0; + TParticle *mPart=0; for (Int_t part = 0; part < npart; part++) { fTrackList[part] = 0; - MPart = gAlice->Particle(part); - mpart = MPart->GetPdgCode(); - child1 = MPart->GetFirstDaughter(); - child2 = MPart->GetLastDaughter(); - pT = MPart->Pt(); - phi = MPart->Phi(); - eta = MPart->Eta(); + mPart = gAlice->GetMCApp()->Particle(part); + mpart = mPart->GetPdgCode(); + child1 = mPart->GetFirstDaughter(); + child2 = mPart->GetLastDaughter(); + pT = mPart->Pt(); + phi = mPart->Phi(); + eta = mPart->Eta(); - px = MPart->Px(); - py = MPart->Py(); - pz = MPart->Pz(); - e = MPart->Energy(); + px = mPart->Px(); + py = mPart->Py(); + pz = mPart->Pz(); + e = mPart->Energy(); // see pyedit in Pythia's text geantPdg = mpart; @@ -1158,10 +1251,10 @@ void AliEMCALJetFinder::FillFromParticles() // printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n", // part, mpart, geantPdg, px, py, pz, e, child1, name.Data()); - PX += px; - PY += py; - PZ += pz; - E += e; + pX += px; + pY += py; + pZ += pz; + energy += e; if (TMath::Abs(eta) <= 0.9) fNChTpc++; @@ -1177,7 +1270,7 @@ void AliEMCALJetFinder::FillFromParticles() } // primary loop printf("\n PX %8.2f PY %8.2f PZ %8.2f E %8.2f \n", - PX, PY, PZ, E); + pX, pY, pZ, energy); DumpLego(); if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc); } @@ -1204,15 +1297,15 @@ void AliEMCALJetFinder::FillFromPartons() // Go through the partons Int_t statusCode=0; for (Int_t part = 8; part < npart; part++) { - TParticle *MPart = gAlice->Particle(part); - Int_t mpart = MPart->GetPdgCode(); + TParticle *mPart = gAlice->GetMCApp()->Particle(part); + Int_t mpart = mPart->GetPdgCode(); // Int_t child1 = MPart->GetFirstDaughter(); - Float_t pT = MPart->Pt(); + Float_t pT = mPart->Pt(); // Float_t p = MPart->P(); - Float_t phi = MPart->Phi(); - Float_t eta = MPart->Eta(); + Float_t phi = mPart->Phi(); + Float_t eta = mPart->Eta(); // Float_t theta = MPart->Theta(); - statusCode = MPart->GetStatusCode(); + statusCode = mPart->GetStatusCode(); // accept partons (21 - gluon, 92 - string) if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue; @@ -1245,8 +1338,8 @@ void AliEMCALJetFinder::BuildTrackFlagTable() { // 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) + TTree *treeK = gAlice->TreeK(); // Get the number of entries in the kine tree + Int_t nKTrks = (Int_t) treeK->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 @@ -1255,8 +1348,10 @@ void AliEMCALJetFinder::BuildTrackFlagTable() { 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 // @@ -1278,7 +1373,7 @@ void AliEMCALJetFinder::BuildTrackFlagTable() { 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; @@ -1295,7 +1390,7 @@ void AliEMCALJetFinder::BuildTrackFlagTable() { //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; @@ -1312,7 +1407,7 @@ void AliEMCALJetFinder::BuildTrackFlagTable() { } //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; @@ -1332,7 +1427,7 @@ void AliEMCALJetFinder::BuildTrackFlagTable() { Int_t AliEMCALJetFinder ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) { - +// Set the flag for the track Int_t flag = 0; Int_t parton = 0; Int_t neutral = 0; @@ -1351,14 +1446,14 @@ Int_t AliEMCALJetFinder -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()); } @@ -1387,7 +1482,20 @@ void AliEMCALJetFinder::SaveBackgroundEvent() 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; iGetSize(); 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() @@ -1400,7 +1508,7 @@ 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"); @@ -1529,14 +1637,11 @@ Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curl // Propagates phi angle to EMCAL radius // static Float_t b = 0.0, rEMCAL = -1.0; - if(rEMCAL<0) { // Get field in kGS - b = gAlice->Field()->SolenoidField(); + b = gAlice->Field()->SolenoidField(); // Get EMCAL radius in cm - rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance(); - printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL); - } - Float_t dPhi = 0.; + rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance(); + Float_t dPhi = 0.; // // // bending radies @@ -1559,20 +1664,33 @@ Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curl return dPhi; } -void hf1(Int_t& id, Float_t& x, Float_t& wgt) +void hf1(Int_t& , Float_t& , Float_t& ) { // dummy for hbook calls ; } void AliEMCALJetFinder::DrawLego(Char_t *opt) -{fLego->Draw(opt);} +{ +// Draw lego + if(fLego) fLego->Draw(opt); +} + +void AliEMCALJetFinder::DrawLegoBackground(Char_t *opt) +{ +// Draw background lego + if(fLegoB) fLegoB->Draw(opt); +} void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt) -{fhLegoEMCAL->Draw(opt);} +{ +// Draw EMCAL Lego + if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt); +} void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode) { +// Draw all hists static TPaveText *varLabel=0; if(!fC1) { fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800); @@ -1613,6 +1731,8 @@ void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode) void AliEMCALJetFinder::PrintParameters(Int_t mode) { +// Print all parameters out + FILE *file=0; if(mode==0) file = stdout; // output to terminal else { @@ -1620,6 +1740,7 @@ void AliEMCALJetFinder::PrintParameters(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); @@ -1637,11 +1758,14 @@ void AliEMCALJetFinder::PrintParameters(Int_t mode) 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); } void AliEMCALJetFinder::DrawLegos() { +// Draw all legos if(!fC1) { fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800); } @@ -1679,6 +1803,7 @@ void AliEMCALJetFinder::DrawLegos() const Char_t* AliEMCALJetFinder::GetFileNameForParameters(Char_t* dir) { +// Get paramters from a file static TString tmp; if(strlen(dir)) tmp = dir; tmp += GetTitle(); @@ -1707,6 +1832,7 @@ void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart) Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg) { +// Return quark info Int_t absPdg = TMath::Abs(pdg); if(absPdg<=6) return kTRUE; // quarks if(pdg == 21) return kTRUE; // gluon @@ -1848,3 +1974,35 @@ void AliEMCALJetFinder::FindChargedJet() if (fWrite) WriteJets(); fEvent++; } +// 16-jan-2003 - just for convenience +void AliEMCALJetFinder::Browse(TBrowser* b) +{ +// Add to browser + if(fHistsList) b->Add((TObject*)fHistsList); +} + +Bool_t AliEMCALJetFinder::IsFolder() const +{ +// Return folder status + 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(" Name of variant %s \n", nt.Data()); + return nt.Data(); +}