From abde65b83b8afb95b417e296477eedef5b9e3c16 Mon Sep 17 00:00:00 2001 From: gconesab Date: Fri, 3 Jul 2009 06:38:17 +0000 Subject: [PATCH] New analysis for electron identification --- PWG4/CMake_libPWG4PartCorrDep.txt | 3 +- PWG4/PWG4PartCorrDepLinkDef.h | 1 + PWG4/PartCorrBase/AliMCAnalysisUtils.cxx | 775 ++++++++++-------- PWG4/PartCorrBase/AliMCAnalysisUtils.h | 120 +-- PWG4/PartCorrDep/AliAnaElectron.cxx | 726 ++++++++++++++++ PWG4/PartCorrDep/AliAnaElectron.h | 177 ++++ PWG4/libPWG4PartCorrDep.pkg | 3 +- .../macros/electrons/ConfigAnalysisElectron.C | 89 ++ PWG4/macros/electrons/GetProofLogs.C | 7 + PWG4/macros/electrons/README | 538 ++++++++++++ PWG4/macros/electrons/anaElectron.C | 535 ++++++++++++ PWG4/macros/electrons/anaElectron.jdl | 34 + PWG4/macros/electrons/anaElectron.sh | 15 + PWG4/macros/electrons/mergeElectron.jdl | 12 + PWG4/macros/electrons/mylauncher.C | 101 +++ PWG4/macros/electrons/parmaker | 104 +++ 16 files changed, 2824 insertions(+), 416 deletions(-) create mode 100755 PWG4/PartCorrDep/AliAnaElectron.cxx create mode 100755 PWG4/PartCorrDep/AliAnaElectron.h create mode 100644 PWG4/macros/electrons/ConfigAnalysisElectron.C create mode 100644 PWG4/macros/electrons/GetProofLogs.C create mode 100644 PWG4/macros/electrons/README create mode 100644 PWG4/macros/electrons/anaElectron.C create mode 100644 PWG4/macros/electrons/anaElectron.jdl create mode 100644 PWG4/macros/electrons/anaElectron.sh create mode 100644 PWG4/macros/electrons/mergeElectron.jdl create mode 100644 PWG4/macros/electrons/mylauncher.C create mode 100755 PWG4/macros/electrons/parmaker diff --git a/PWG4/CMake_libPWG4PartCorrDep.txt b/PWG4/CMake_libPWG4PartCorrDep.txt index 88722b3b84c..68eea3a9618 100755 --- a/PWG4/CMake_libPWG4PartCorrDep.txt +++ b/PWG4/CMake_libPWG4PartCorrDep.txt @@ -7,7 +7,8 @@ set(SRCS PartCorrDep/AliAnaParticleIsolation.cxx PartCorrDep/AliAnaParticlePartonCorrelation.cxx PartCorrDep/AliAnaParticleHadronCorrelation.cxx PartCorrDep/AliAnaParticleJetFinderCorrelation.cxx PartCorrDep/AliAnaParticleJetLeadingConeCorrelation.cxx PartCorrDep/AliAnaChargedParticles.cxx - PartCorrDep/AliAnaCalorimeterQA.cxx PartCorrDep/AliAnaNeutralMeson.cxx + PartCorrDep/AliAnaCalorimeterQA.cxx PartCorrDep/AliAnaNeutralMeson.cxx + PartCorrDep/AliAnaElectron.cxx ) # fill list of header files from list of source files diff --git a/PWG4/PWG4PartCorrDepLinkDef.h b/PWG4/PWG4PartCorrDepLinkDef.h index 27d79c49ac9..9c2737dfdaa 100755 --- a/PWG4/PWG4PartCorrDepLinkDef.h +++ b/PWG4/PWG4PartCorrDepLinkDef.h @@ -18,5 +18,6 @@ #pragma link C++ class AliAnaParticleJetLeadingConeCorrelation+; #pragma link C++ class AliAnaCalorimeterQA+; #pragma link C++ class AliAnaNeutralMeson+; +#pragma link C++ class AliAnaElectron+; #endif diff --git a/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx b/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx index c02dbfd9bd9..b1a567e2a9e 100755 --- a/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx +++ b/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx @@ -1,354 +1,421 @@ -/************************************************************************** - * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * * - * Author: The ALICE Off-line Project. * - * Contributors are mentioned in the code where appropriate. * - * * - * Permission to use, copy, modify and distribute this software and its * - * documentation strictly for non-commercial purposes is hereby granted * - * without fee, provided that the above copyright notice appears in all * - * copies and that both the copyright notice and this permission notice * - * appear in the supporting documentation. The authors make no claims * - * about the suitability of this software for any purpose. It is * - * provided "as is" without express or implied warranty. * - **************************************************************************/ -/* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */ - -//_________________________________________________________________________ -// Class for analysis utils for MC data -// stored in stack or event header. -// Contains: -// - method to check the origin of a given track/cluster -// - method to obtain the generated jets -// -//*-- Author: Gustavo Conesa (LNF-INFN) -////////////////////////////////////////////////////////////////////////////// - - -// --- ROOT system --- -#include -#include - -//---- ANALYSIS system ---- -#include "AliMCAnalysisUtils.h" -#include "AliStack.h" -#include "TParticle.h" -#include "AliGenPythiaEventHeader.h" - - ClassImp(AliMCAnalysisUtils) - - //________________________________________________ - AliMCAnalysisUtils::AliMCAnalysisUtils() : - TObject(), fCurrentEvent(-1), fDebug(-1), - fJetsList(new TList), fMCGenerator("PYTHIA") -{ - //Ctor -} - -//____________________________________________________________________________ -AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) : - TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug), - fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator) -{ - // cpy ctor - -} - -//_________________________________________________________________________ -AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils) -{ - // assignment operator - - if(&mcutils == this) return *this; - fCurrentEvent = mcutils.fCurrentEvent ; - fDebug = mcutils.fDebug; - fJetsList = mcutils.fJetsList; - fMCGenerator = mcutils.fMCGenerator; - - return *this; -} - -//____________________________________________________________________________ -AliMCAnalysisUtils::~AliMCAnalysisUtils() -{ - // Remove all pointers. - - if (fJetsList) { - fJetsList->Clear(); - delete fJetsList ; - } -} - -//_________________________________________________________________________ -Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliStack * stack) const { - //Play with the MC stack if available - //Check origin of the candidates, good for PYTHIA - - if(!stack) { - printf("AliMCAnalysisUtils::CheckOrigin() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); - abort(); - } - // printf("label %d, ntrack %d, nprim %d\n",label, stack->GetNtrack(), stack->GetNprimary()); - // for(Int_t i = 0; i< stack->GetNprimary(); i++){ - // TParticle *particle = stack->Particle(i); - // //particle->Print(); - // } - if(label >= 0 && label < stack->GetNtrack()){ - //Mother - TParticle * mom = stack->Particle(label); - Int_t mPdg = TMath::Abs(mom->GetPdgCode()); - Int_t mStatus = mom->GetStatusCode() ; - Int_t iParent = mom->GetFirstMother() ; - if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOrigin: Mother is parton %d\n",iParent); - - //GrandParent - TParticle * parent = new TParticle ; - Int_t pPdg = -1; - Int_t pStatus =-1; - if(iParent > 0){ - parent = stack->Particle(iParent); - pPdg = TMath::Abs(parent->GetPdgCode()); - pStatus = parent->GetStatusCode(); - } - else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOrigin: Parent with label %d\n",iParent); - - //return tag - if(mPdg == 22){ - if(mStatus == 1){ - if(fMCGenerator == "PYTHIA"){ - if(iParent < 8 && iParent > 5) {//outgoing partons - if(pPdg == 22) return kMCPrompt; - else return kMCFragmentation; - }//Outgoing partons - else if(pStatus == 11){//Decay - if(pPdg == 111) return kMCPi0Decay ; - else if (pPdg == 221) return kMCEtaDecay ; - else return kMCOtherDecay ; - }//Decay - else return kMCISR; //Initial state radiation - }//PYTHIA - - else if(fMCGenerator == "HERWIG"){ - if(pStatus < 197){//Not decay - while(1){ - if(parent->GetFirstMother()<=5) break; - iParent = parent->GetFirstMother(); - parent=stack->Particle(iParent); - pStatus= parent->GetStatusCode(); - pPdg = parent->GetPdgCode(); - }//Look for the parton - - if(iParent < 8 && iParent > 5) { - if(pPdg == 22) return kMCPrompt; - else return kMCFragmentation; - } - return kMCISR;//Initial state radiation - }//Not decay - else{//Decay - if(pPdg == 111) return kMCPi0Decay ; - else if (pPdg == 221) return kMCEtaDecay ; - else return kMCOtherDecay ; - }//Decay - }//HERWIG - else return kMCUnknown; - }//Status 1 : Pythia generated - else if(mStatus == 0){ - if(pPdg ==22 || pPdg ==11|| pPdg == 2112 || pPdg == 211 || - pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 ) - return kMCConversion ; - if(pPdg == 111) return kMCPi0Decay ; - else if (pPdg == 221) return kMCEtaDecay ; - else return kMCOtherDecay ; - }//status 0 : geant generated - }//Mother Photon - else if(mPdg == 111) return kMCPi0 ; - else if(mPdg == 221) return kMCEta ; - else if(mPdg ==11){ -// if(pPdg !=22 &&mStatus == 0) { -// printf("Origin electron, pT %f, status %d, parent %s, pstatus %d, vx %f, vy %f, vz %f\n", -// mom->Pt(),mStatus, parent->GetName(),pStatus,mom->Vx(),mom->Vy(), mom->Vz()); - -// } - - if(mStatus == 0) { - if(pPdg ==22) return kMCConversion ; - else if (pStatus == 0) return kMCConversion; - else return kMCElectron ; - } - else return kMCElectron ; - } - else return kMCUnknown; - }//Good label value - else{ - if(label < 0 ) printf("AliMCAnalysisUtils::CheckOrigin: *** bad label or no stack ***: label %d \n", label); - if(label >= stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOrigin: *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); - return kMCUnknown; - }//Bad label - - return kMCUnknown; - -} - -//_________________________________________________________________________ -TList * AliMCAnalysisUtils::GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) { - //Return list of jets (TParticles) and index of most likely parton that originated it. - - if(fCurrentEvent!=iEvent){ - fCurrentEvent = iEvent; - fJetsList = new TList; - Int_t nTriggerJets = 0; - Float_t tmpjet[]={0,0,0,0}; - - //printf("Event %d %d\n",fCurrentEvent,iEvent); - //Get outgoing partons - if(stack->GetNtrack() < 8) return fJetsList; - TParticle * parton1 = stack->Particle(6); - TParticle * parton2 = stack->Particle(7); - if(fDebug > 2){ - printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", - parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()); - printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", - parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()); - } -// //Trace the jet from the mother parton -// Float_t pt = 0; -// Float_t pt1 = 0; -// Float_t pt2 = 0; -// Float_t e = 0; -// Float_t e1 = 0; -// Float_t e2 = 0; -// TParticle * tmptmp = new TParticle; -// for(Int_t i = 0; i< stack->GetNprimary(); i++){ -// tmptmp = stack->Particle(i); - -// if(tmptmp->GetStatusCode() == 1){ -// pt = tmptmp->Pt(); -// e = tmptmp->Energy(); -// Int_t imom = tmptmp->GetFirstMother(); -// Int_t imom1 = 0; -// //printf("1st imom %d\n",imom); -// while(imom > 5){ -// imom1=imom; -// tmptmp = stack->Particle(imom); -// imom = tmptmp->GetFirstMother(); -// //printf("imom %d \n",imom); -// } -// //printf("Last imom %d %d\n",imom1, imom); -// if(imom1 == 6) { -// pt1+=pt; -// e1+=e; -// } -// else if (imom1 == 7){ -// pt2+=pt; -// e2+=e; } -// }// status 1 - -// }// for - -// printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2); - - //Get the jet, different way for different generator - //PYTHIA - if(fMCGenerator == "PYTHIA"){ - TParticle * jet = new TParticle; - AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh; - nTriggerJets = pygeh->NTriggerJets(); - if(fDebug > 1) - printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets); - - Int_t iparton = -1; - for(Int_t i = 0; i< nTriggerJets; i++){ - iparton=-1; - pygeh->TriggerJet(i, tmpjet); - jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0); - //Assign an outgoing parton as mother - Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi()); - Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi()); - if(phidiff1 > phidiff2) jet->SetFirstMother(7); - else jet->SetFirstMother(6); - //jet->Print(); - if(fDebug > 1) - printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", - i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()); - fJetsList->Add(jet); - } - }//Pythia triggered jets - //HERWIG - else if (fMCGenerator=="HERWIG"){ - Int_t pdg = -1; - //Check parton 1 - TParticle * tmp = parton1; - if(parton1->GetPdgCode()!=22){ - while(pdg != 94){ - if(tmp->GetFirstDaughter()==-1) return fJetsList; - tmp = stack->Particle(tmp->GetFirstDaughter()); - pdg = tmp->GetPdgCode(); - }//while - - //Add found jet to list - TParticle *jet1 = new TParticle(*tmp); - jet1->SetFirstMother(6); - fJetsList->Add(jet1); - //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter()); - //tmp = stack->Particle(tmp->GetFirstDaughter()); - //tmp->Print(); - //jet1->Print(); - if(fDebug > 1) - printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", - jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()); - }//not photon - - //Check parton 2 - pdg = -1; - tmp = parton2; - Int_t i = -1; - if(parton2->GetPdgCode()!=22){ - while(pdg != 94){ - if(tmp->GetFirstDaughter()==-1) return fJetsList; - i = tmp->GetFirstDaughter(); - tmp = stack->Particle(tmp->GetFirstDaughter()); - pdg = tmp->GetPdgCode(); - }//while - //Add found jet to list - TParticle *jet2 = new TParticle(*tmp); - jet2->SetFirstMother(7); - fJetsList->Add(jet2); - //jet2->Print(); - if(fDebug > 1) - printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", - jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()); - //Int_t first = tmp->GetFirstDaughter(); - //Int_t last = tmp->GetLastDaughter(); - //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode()); - // for(Int_t d = first ; d < last+1; d++){ -// tmp = stack->Particle(d); -// if(i == tmp->GetFirstMother()) -// printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", -// d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta()); -// } - //tmp->Print(); - }//not photon - }//Herwig generated jets - } - - return fJetsList; -} - -//________________________________________________________________ -void AliMCAnalysisUtils::Print(const Option_t * opt) const -{ - //Print some relevant parameters set for the analysis - - if(! opt) - return; - - printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ; - - printf("Debug level = %d\n",fDebug); - printf("MC Generator = %s\n",fMCGenerator.Data()); - - printf(" \n"); - -} - - +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ +/* $Id: AliMCAnalysisUtils.cxx 21839 2007-10-29 13:49:42Z gustavo $ */ + +//_________________________________________________________________________ +// Class for analysis utils for MC data +// stored in stack or event header. +// Contains: +// - method to check the origin of a given track/cluster +// - method to obtain the generated jets +// +//*-- Author: Gustavo Conesa (LNF-INFN) +////////////////////////////////////////////////////////////////////////////// + + +// --- ROOT system --- +#include +#include + +//---- ANALYSIS system ---- +#include "AliMCAnalysisUtils.h" +#include "AliStack.h" +#include "TParticle.h" +#include "AliGenPythiaEventHeader.h" + + ClassImp(AliMCAnalysisUtils) + + //________________________________________________ + AliMCAnalysisUtils::AliMCAnalysisUtils() : + TObject(), fCurrentEvent(-1), fDebug(-1), + fJetsList(new TList), fMCGenerator("PYTHIA") +{ + //Ctor +} + +//____________________________________________________________________________ +AliMCAnalysisUtils::AliMCAnalysisUtils(const AliMCAnalysisUtils & mcutils) : + TObject(mcutils), fCurrentEvent(mcutils.fCurrentEvent), fDebug(mcutils.fDebug), + fJetsList(mcutils.fJetsList), fMCGenerator(mcutils.fMCGenerator) +{ + // cpy ctor + +} + +//_________________________________________________________________________ +AliMCAnalysisUtils & AliMCAnalysisUtils::operator = (const AliMCAnalysisUtils & mcutils) +{ + // assignment operator + + if(&mcutils == this) return *this; + fCurrentEvent = mcutils.fCurrentEvent ; + fDebug = mcutils.fDebug; + fJetsList = mcutils.fJetsList; + fMCGenerator = mcutils.fMCGenerator; + + return *this; +} + +//____________________________________________________________________________ +AliMCAnalysisUtils::~AliMCAnalysisUtils() +{ + // Remove all pointers. + + if (fJetsList) { + fJetsList->Clear(); + delete fJetsList ; + } +} + +//_________________________________________________________________________ +Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliStack * stack) const { + //Play with the MC stack if available + //Check origin of the candidates, good for PYTHIA + + if(!stack) { + printf("AliMCAnalysisUtils::CheckOrigin() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); + abort(); + } + // printf("label %d, ntrack %d, nprim %d\n",label, stack->GetNtrack(), stack->GetNprimary()); + // for(Int_t i = 0; i< stack->GetNprimary(); i++){ + // TParticle *particle = stack->Particle(i); + // //particle->Print(); + // } + if(label >= 0 && label < stack->GetNtrack()){ + //Mother + TParticle * mom = stack->Particle(label); + Int_t mPdg = TMath::Abs(mom->GetPdgCode()); + Int_t mStatus = mom->GetStatusCode() ; + Int_t iParent = mom->GetFirstMother() ; + if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOrigin: Mother is parton %d\n",iParent); + + //GrandParent + TParticle * parent = new TParticle ; + Int_t pPdg = -1; + Int_t pStatus =-1; + if(iParent > 0){ + parent = stack->Particle(iParent); + pPdg = TMath::Abs(parent->GetPdgCode()); + pStatus = parent->GetStatusCode(); + } + else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOrigin: Parent with label %d\n",iParent); + + //return tag + if(mPdg == 22){ //photon + if(mStatus == 1){ //undecayed particle + if(fMCGenerator == "PYTHIA"){ + if(iParent < 8 && iParent > 5) {//outgoing partons + if(pPdg == 22) return kMCPrompt; + else return kMCFragmentation; + }//Outgoing partons + else if(pStatus == 11){//Decay + if(pPdg == 111) return kMCPi0Decay ; + else if (pPdg == 221) return kMCEtaDecay ; + else return kMCOtherDecay ; + }//Decay + else return kMCISR; //Initial state radiation + }//PYTHIA + + else if(fMCGenerator == "HERWIG"){ + if(pStatus < 197){//Not decay + while(1){ + if(parent->GetFirstMother()<=5) break; + iParent = parent->GetFirstMother(); + parent=stack->Particle(iParent); + pStatus= parent->GetStatusCode(); + pPdg = parent->GetPdgCode(); + }//Look for the parton + + if(iParent < 8 && iParent > 5) { + if(pPdg == 22) return kMCPrompt; + else return kMCFragmentation; + } + return kMCISR;//Initial state radiation + }//Not decay + else{//Decay + if(pPdg == 111) return kMCPi0Decay ; + else if (pPdg == 221) return kMCEtaDecay ; + else return kMCOtherDecay ; + }//Decay + }//HERWIG + else return kMCUnknown; + }//Status 1 : Pythia generated + else if(mStatus == 0){ + if(pPdg ==22 || pPdg ==11|| pPdg == 2112 || pPdg == 211 || + pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 ) + return kMCConversion ; + if(pPdg == 111) return kMCPi0Decay ; + else if (pPdg == 221) return kMCEtaDecay ; + else return kMCOtherDecay ; + }//status 0 : geant generated + }//Mother Photon + else if(mPdg == 111) return kMCPi0 ; + else if(mPdg == 221) return kMCEta ; + + //cluster's mother is an electron. Where did that electron come from? + else if(mPdg == 11){ //electron + + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOrigin: Checking ancestors of electrons"); + + //check first for B and C ancestry, then other possibilities. + //An electron from a photon parent could have other particles in + //its history and we would want to know that, right? + + if(mStatus == 1) { //electron from event generator + if (pPdg == -1) return kMCElectron; //no parent + else if (pPdg == 23) return kMCZDecay; //parent is Z-boson + else if (pPdg == 24) return kMCWDecay; //parent is W-boson + else { //check the electron's ancestors for B/C contribution + Bool_t BAncestor = kFALSE; + Bool_t CAncestor = kFALSE; + TParticle * ancestors = stack->Particle(label); + Int_t aPdg = TMath::Abs(ancestors->GetPdgCode()); + //Int_t aStatus = ancestors->GetStatusCode(); + Int_t iAncestors = ancestors->GetFirstMother(); + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOrigin: Scaning the decay chain for bottom/charm generated electron"); + while(ancestors->IsPrimary()){//searching for ancestors + if((499 < aPdg && aPdg < 600)||(4999 < aPdg && aPdg < 6000)) BAncestor = kTRUE; + if((399 < aPdg && aPdg < 500)||(3999 < aPdg && aPdg < 5000)) CAncestor = kTRUE; + if(BAncestor && CAncestor) break; + iAncestors = ancestors->GetFirstMother(); + ancestors = stack->Particle(iAncestors); + aPdg = ancestors->GetPdgCode(); + }//searching for ancestors + if(BAncestor && CAncestor) return kMCEFromCFromB;//Decay chain has both B and C + else if(BAncestor && !CAncestor) return kMCEFromB;//Decay chain has only B + else if(!BAncestor && CAncestor) return kMCEFromC;//Decay chain has only C + } + //if it is not from W,Z or B/C ancestor, where is it from? + if (pPdg == 111) return kMCPi0Decay;//Pi0 Dalitz decay + else if(pPdg == 221) return kMCEtaDecay;//Eta Dalitz decay + else return kMCOtherDecay; + + } else if (mStatus == 0) { //electron from GEANT + + //Rewind ancestry and check for electron with status == 1 + //if we find one, we'll assume that this object is from an + //electron but that it may have gone through some showering in + //material before the detector + + //Not a double-counting problem because we are only accessing + //these histories for MC labels connected to a reco object. + //If you wanted to use this to sort through the kine stack + //directly, might it be a problem? + Bool_t EleFromEvGen = kFALSE; + Bool_t BAncestor = kFALSE; + Bool_t CAncestor = kFALSE; + + TParticle * ancestors = stack->Particle(label); + Int_t aPdg = TMath::Abs(ancestors->GetPdgCode()); + Int_t aStatus = ancestors->GetStatusCode(); + Int_t iAncestors = ancestors->GetFirstMother(); + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOrigin: Scaning the decay chain for bottom/charm electrons"); + while(ancestors->IsPrimary()){//searching for ancestors + if(aStatus == 1 && aPdg == 11) EleFromEvGen = kTRUE; + if(EleFromEvGen && aPdg == 23) return kMCZDecay; + if(EleFromEvGen && aPdg == 24) return kMCWDecay; + if(EleFromEvGen && ((499 < aPdg && aPdg < 600)||(4999 < aPdg && aPdg < 6000))) BAncestor = kTRUE; + if(EleFromEvGen && ((399 < aPdg && aPdg < 500)||(3999 < aPdg && aPdg < 5000))) CAncestor = kTRUE; + if(BAncestor && CAncestor) break; + iAncestors = ancestors->GetFirstMother(); + ancestors = stack->Particle(iAncestors); + aPdg = ancestors->GetPdgCode(); + }//searching for ancestors + if(BAncestor && CAncestor) return kMCEFromCFromB;//Decay chain has both B and C + else if(BAncestor && !CAncestor) return kMCEFromB;//Decay chain has only B + else if(!BAncestor && CAncestor) return kMCEFromC;//Decay chain has only C + if(pPdg ==22 || pPdg ==11|| pPdg == 2112 || pPdg == 211 || + pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 ) + return kMCConversion ; + if(pPdg == 111) return kMCPi0Decay ; + else if (pPdg == 221) return kMCEtaDecay ; + else return kMCOtherDecay ; + } //GEANT check + }//electron check + else return kMCUnknown; + }//Good label value + else{ + if(label < 0 ) printf("AliMCAnalysisUtils::CheckOrigin: *** bad label or no stack ***: label %d \n", label); + if(label >= stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOrigin: *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); + return kMCUnknown; + }//Bad label + + return kMCUnknown; + +} + +//_________________________________________________________________________ +TList * AliMCAnalysisUtils::GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) { + //Return list of jets (TParticles) and index of most likely parton that originated it. + + if(fCurrentEvent!=iEvent){ + fCurrentEvent = iEvent; + fJetsList = new TList; + Int_t nTriggerJets = 0; + Float_t tmpjet[]={0,0,0,0}; + + //printf("Event %d %d\n",fCurrentEvent,iEvent); + //Get outgoing partons + if(stack->GetNtrack() < 8) return fJetsList; + TParticle * parton1 = stack->Particle(6); + TParticle * parton2 = stack->Particle(7); + if(fDebug > 2){ + printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", + parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()); + printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", + parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()); + } +// //Trace the jet from the mother parton +// Float_t pt = 0; +// Float_t pt1 = 0; +// Float_t pt2 = 0; +// Float_t e = 0; +// Float_t e1 = 0; +// Float_t e2 = 0; +// TParticle * tmptmp = new TParticle; +// for(Int_t i = 0; i< stack->GetNprimary(); i++){ +// tmptmp = stack->Particle(i); + +// if(tmptmp->GetStatusCode() == 1){ +// pt = tmptmp->Pt(); +// e = tmptmp->Energy(); +// Int_t imom = tmptmp->GetFirstMother(); +// Int_t imom1 = 0; +// //printf("1st imom %d\n",imom); +// while(imom > 5){ +// imom1=imom; +// tmptmp = stack->Particle(imom); +// imom = tmptmp->GetFirstMother(); +// //printf("imom %d \n",imom); +// } +// //printf("Last imom %d %d\n",imom1, imom); +// if(imom1 == 6) { +// pt1+=pt; +// e1+=e; +// } +// else if (imom1 == 7){ +// pt2+=pt; +// e2+=e; } +// }// status 1 + +// }// for + +// printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2); + + //Get the jet, different way for different generator + //PYTHIA + if(fMCGenerator == "PYTHIA"){ + TParticle * jet = new TParticle; + AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh; + nTriggerJets = pygeh->NTriggerJets(); + if(fDebug > 1) + printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets); + + Int_t iparton = -1; + for(Int_t i = 0; i< nTriggerJets; i++){ + iparton=-1; + pygeh->TriggerJet(i, tmpjet); + jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0); + //Assign an outgoing parton as mother + Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi()); + Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi()); + if(phidiff1 > phidiff2) jet->SetFirstMother(7); + else jet->SetFirstMother(6); + //jet->Print(); + if(fDebug > 1) + printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", + i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()); + fJetsList->Add(jet); + } + }//Pythia triggered jets + //HERWIG + else if (fMCGenerator=="HERWIG"){ + Int_t pdg = -1; + //Check parton 1 + TParticle * tmp = parton1; + if(parton1->GetPdgCode()!=22){ + while(pdg != 94){ + if(tmp->GetFirstDaughter()==-1) return fJetsList; + tmp = stack->Particle(tmp->GetFirstDaughter()); + pdg = tmp->GetPdgCode(); + }//while + + //Add found jet to list + TParticle *jet1 = new TParticle(*tmp); + jet1->SetFirstMother(6); + fJetsList->Add(jet1); + //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter()); + //tmp = stack->Particle(tmp->GetFirstDaughter()); + //tmp->Print(); + //jet1->Print(); + if(fDebug > 1) + printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", + jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()); + }//not photon + + //Check parton 2 + pdg = -1; + tmp = parton2; + Int_t i = -1; + if(parton2->GetPdgCode()!=22){ + while(pdg != 94){ + if(tmp->GetFirstDaughter()==-1) return fJetsList; + i = tmp->GetFirstDaughter(); + tmp = stack->Particle(tmp->GetFirstDaughter()); + pdg = tmp->GetPdgCode(); + }//while + //Add found jet to list + TParticle *jet2 = new TParticle(*tmp); + jet2->SetFirstMother(7); + fJetsList->Add(jet2); + //jet2->Print(); + if(fDebug > 1) + printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", + jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()); + //Int_t first = tmp->GetFirstDaughter(); + //Int_t last = tmp->GetLastDaughter(); + //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode()); + // for(Int_t d = first ; d < last+1; d++){ +// tmp = stack->Particle(d); +// if(i == tmp->GetFirstMother()) +// printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", +// d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta()); +// } + //tmp->Print(); + }//not photon + }//Herwig generated jets + } + + return fJetsList; +} + +//________________________________________________________________ +void AliMCAnalysisUtils::Print(const Option_t * opt) const +{ + //Print some relevant parameters set for the analysis + + if(! opt) + return; + + printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ; + + printf("Debug level = %d\n",fDebug); + printf("MC Generator = %s\n",fMCGenerator.Data()); + + printf(" \n"); + +} + + diff --git a/PWG4/PartCorrBase/AliMCAnalysisUtils.h b/PWG4/PartCorrBase/AliMCAnalysisUtils.h index 22f0bbbd653..a0f5fc1d79e 100755 --- a/PWG4/PartCorrBase/AliMCAnalysisUtils.h +++ b/PWG4/PartCorrBase/AliMCAnalysisUtils.h @@ -1,60 +1,60 @@ -#ifndef ALIMCANALYSISUTILS_H -#define ALIMCANALYSISUTILS_H -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ -/* $Id: $ */ - -//_________________________________________________________________________ -// Class for analysis utils for MC data -// stored in stack or event header. -// Contains: -// - method to check the origin of a given track/cluster -// - method to obtain the generated jets -// -//*-- Author: Gustavo Conesa (INFN-LNF) - -// --- ROOT system --- -#include -class TString ; -class TList ; - -//--- AliRoot system --- -class AliStack ; -class AliGenEventHeader ; - -class AliMCAnalysisUtils : public TObject { - -public: - - AliMCAnalysisUtils() ; // ctor - AliMCAnalysisUtils(const AliMCAnalysisUtils & g) ; // cpy ctor - AliMCAnalysisUtils & operator = (const AliMCAnalysisUtils & g) ;//cpy assignment - virtual ~AliMCAnalysisUtils() ;//virtual dtor - - enum mcTypes {kMCPrompt, kMCFragmentation, kMCISR, kMCPi0Decay, kMCEtaDecay, kMCOtherDecay, kMCPi0, kMCEta, kMCElectron, kMCConversion, kMCUnknown}; - - Int_t CheckOrigin(const Int_t label, AliStack * stack) const ; - TList * GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) ; - - void Print(const Option_t * opt)const; - - void SetDebug(Int_t deb) {fDebug=deb;} - Int_t GetDebug() const {return fDebug;} - - void SetMCGenerator(TString mcgen) {fMCGenerator=mcgen;} - TString GetMCGenerator() const {return fMCGenerator;} - -private: - Int_t fCurrentEvent; // Current Event - Int_t fDebug; // Debug level - TList * fJetsList; // List of jets - TString fMCGenerator; // MC geneator used to generate data in simulation - - ClassDef(AliMCAnalysisUtils,1) -} ; - - -#endif //ALIMCANALYSISUTILS_H - - - +#ifndef ALIMCANALYSISUTILS_H +#define ALIMCANALYSISUTILS_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +/* $Id: $ */ + +//_________________________________________________________________________ +// Class for analysis utils for MC data +// stored in stack or event header. +// Contains: +// - method to check the origin of a given track/cluster +// - method to obtain the generated jets +// +//*-- Author: Gustavo Conesa (INFN-LNF) + +// --- ROOT system --- +#include +class TString ; +class TList ; + +//--- AliRoot system --- +class AliStack ; +class AliGenEventHeader ; + +class AliMCAnalysisUtils : public TObject { + +public: + + AliMCAnalysisUtils() ; // ctor + AliMCAnalysisUtils(const AliMCAnalysisUtils & g) ; // cpy ctor + AliMCAnalysisUtils & operator = (const AliMCAnalysisUtils & g) ;//cpy assignment + virtual ~AliMCAnalysisUtils() ;//virtual dtor + + enum mcTypes {kMCPrompt, kMCFragmentation, kMCISR, kMCPi0Decay, kMCEtaDecay, kMCOtherDecay, kMCPi0, kMCEta, kMCElectron, kMCConversion, kMCUnknown, kMCEFromCFromB, kMCEFromC, kMCEFromB,kMCZDecay,kMCWDecay}; + + Int_t CheckOrigin(const Int_t label, AliStack * stack) const ; + TList * GetJets(Int_t iEvent, AliStack * stack, AliGenEventHeader * geh) ; + + void Print(const Option_t * opt)const; + + void SetDebug(Int_t deb) {fDebug=deb;} + Int_t GetDebug() const {return fDebug;} + + void SetMCGenerator(TString mcgen) {fMCGenerator=mcgen;} + TString GetMCGenerator() const {return fMCGenerator;} + +private: + Int_t fCurrentEvent; // Current Event + Int_t fDebug; // Debug level + TList * fJetsList; // List of jets + TString fMCGenerator; // MC geneator used to generate data in simulation + + ClassDef(AliMCAnalysisUtils,1) +} ; + + +#endif //ALIMCANALYSISUTILS_H + + + diff --git a/PWG4/PartCorrDep/AliAnaElectron.cxx b/PWG4/PartCorrDep/AliAnaElectron.cxx new file mode 100755 index 00000000000..6e00ca2fa4c --- /dev/null +++ b/PWG4/PartCorrDep/AliAnaElectron.cxx @@ -0,0 +1,726 @@ + /************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ +/* $Id: $ */ + +//_________________________________________________________________________ +// +// Class for the electron identification. +// Clusters from EMCAL matched to tracks +// and kept in the AOD. Few histograms produced. +// +// -- Author: J.L. Klay (Cal Poly) +////////////////////////////////////////////////////////////////////////////// + + +// --- ROOT system --- +#include +#include +#include +#include + +// --- Analysis system --- +#include "AliAnaElectron.h" +#include "AliCaloTrackReader.h" +#include "AliMCAnalysisUtils.h" +#include "AliFidutialCut.h" +#include "AliESDCaloCluster.h" +//#include "AliESDCaloCells.h" +#include "AliESDtrack.h" +#include "AliESDEvent.h" +#include "AliCaloPID.h" +#include "AliVEvent.h" + +ClassImp(AliAnaElectron) + +//____________________________________________________________________________ +AliAnaElectron::AliAnaElectron() + : AliAnaPartCorrBaseClass(), fCalorimeter(""), + fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.), + //matching checks + fh1pOverE(0),fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0), + fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0), + fh2OuterPtVsExtrapPt(0),fh2OuterPhiVsExtrapPhi(0),fh2OuterEtaVsExtrapEta(0), + fh2TrackPVsClusterE(0),fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0), + //reco + fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0), + //MC + fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0), + fhPtBottom(0),fhPhiBottom(0),fhEtaBottom(0), + fhPtCharm(0),fhPhiCharm(0),fhEtaCharm(0), + fhPtCFromB(0),fhPhiCFromB(0),fhEtaCFromB(0), + fhPtDalitz(0),fhPhiDalitz(0),fhEtaDalitz(0), + fhPtWDecay(0),fhPhiWDecay(0),fhEtaWDecay(0), + fhPtZDecay(0),fhPhiZDecay(0),fhEtaZDecay(0), + fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0), + fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0) +// fhMCElePt(0),fhMCElePhi(0),fhMCEleEta(0) +{ + //default ctor + + //Initialize parameters + InitParameters(); + +} + +//____________________________________________________________________________ +AliAnaElectron::AliAnaElectron(const AliAnaElectron & g) + : AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter), + fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),fResidualCut(g.fResidualCut), + //matching checks + fh1pOverE(g.fh1pOverE),fh1dR(g.fh1dR), + fh2EledEdx(g.fh2EledEdx),fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi), + fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched), + fh2OuterPtVsExtrapPt(g.fh2OuterPtVsExtrapPt),fh2OuterPhiVsExtrapPhi(g.fh2OuterPhiVsExtrapPhi), + fh2OuterEtaVsExtrapEta(g.fh2OuterEtaVsExtrapEta), + fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE), + fh2TrackPhiVsClusterPhi(g.fh2TrackPhiVsClusterPhi),fh2TrackEtaVsClusterEta(g.fh2TrackEtaVsClusterEta), + //reco + fhPtElectron(g.fhPtElectron),fhPhiElectron(g.fhPhiElectron),fhEtaElectron(g.fhEtaElectron), + //MC + fhPtConversion(g.fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion), + fhPtBottom(g.fhPtBottom),fhPhiBottom(g.fhPhiBottom),fhEtaBottom(g.fhEtaBottom), + fhPtCharm(g.fhPtCharm),fhPhiCharm(g.fhPhiCharm),fhEtaCharm(g.fhEtaCharm), + fhPtCFromB(g.fhPtCFromB),fhPhiCFromB(g.fhPhiCFromB),fhEtaCFromB(g.fhEtaCFromB), + fhPtDalitz(g.fhPtDalitz),fhPhiDalitz(g.fhPhiDalitz),fhEtaDalitz(g.fhEtaDalitz), + fhPtWDecay(g.fhPtWDecay),fhPhiWDecay(g.fhPhiWDecay),fhEtaWDecay(g.fhEtaWDecay), + fhPtZDecay(g.fhPtZDecay),fhPhiZDecay(g.fhPhiZDecay),fhEtaZDecay(g.fhEtaZDecay), + fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt), + fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown) + // fhMCElePt(g.fhMCElePt),fhMCElePhi(g.fhMCElePhi),fhMCEleEta(g.fhMCEleEta) +{ + // cpy ctor + +} + +//_________________________________________________________________________ +AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g) +{ + // assignment operator + + if(&g == this) return *this; + fCalorimeter = g.fCalorimeter; + fpOverEmin = g.fpOverEmin; + fpOverEmax = g.fpOverEmax; + fResidualCut = g.fResidualCut; + //fhEnergy = g.fhEnergy; + //fhClusMult = g.fhClusMult; + //fhClusters = g.fhClusters; + //fhDigitsEvent = g.fhDigitsEvent; + fh1pOverE = g.fh1pOverE; + fh1dR = g.fh1dR; + fh2EledEdx = g.fh2EledEdx; + fh2MatchdEdx = g.fh2MatchdEdx; + fh2dEtadPhi = g.fh2dEtadPhi; + fh2dEtadPhiMatched = g.fh2dEtadPhiMatched; + fh2dEtadPhiUnmatched = g.fh2dEtadPhiUnmatched; + fh2OuterPtVsExtrapPt = g.fh2OuterPtVsExtrapPt; + fh2OuterPhiVsExtrapPhi = g.fh2OuterPhiVsExtrapPhi; + fh2OuterEtaVsExtrapEta = g.fh2OuterEtaVsExtrapEta; + fh2TrackPVsClusterE = g.fh2TrackPVsClusterE; + fh2TrackPtVsClusterE = g.fh2TrackPtVsClusterE; + fh2TrackPhiVsClusterPhi = g.fh2TrackPhiVsClusterPhi; + fh2TrackEtaVsClusterEta = g.fh2TrackEtaVsClusterEta; + fhPtElectron = g.fhPtElectron; + fhPhiElectron = g.fhPhiElectron; + fhEtaElectron = g.fhEtaElectron; + fhPtConversion = g.fhPtConversion; + fhPhiConversion = g.fhPhiConversion; + fhEtaConversion = g.fhEtaConversion; + fhPtBottom = g.fhPtBottom; + fhPhiBottom = g.fhPhiBottom; + fhEtaBottom = g.fhEtaBottom; + fhPtCharm = g.fhPtCharm; + fhPhiCharm = g.fhPhiCharm; + fhEtaCharm = g.fhEtaCharm; + fhPtCFromB = g.fhPtCFromB; + fhPhiCFromB = g.fhPhiCFromB; + fhEtaCFromB = g.fhEtaCFromB; + fhPtDalitz = g.fhPtDalitz; + fhPhiDalitz = g.fhPhiDalitz; + fhEtaDalitz = g.fhEtaDalitz; + fhPtWDecay = g.fhPtWDecay; + fhPhiWDecay = g.fhPhiWDecay; + fhEtaWDecay = g.fhEtaWDecay; + fhPtZDecay = g.fhPtZDecay; + fhPhiZDecay = g.fhPhiZDecay; + fhEtaZDecay = g.fhEtaZDecay; + fhPtPrompt = g.fhPtPrompt; + fhPhiPrompt = g.fhPhiPrompt; + fhEtaPrompt = g.fhEtaPrompt; + fhPtUnknown = g.fhPtUnknown; + fhPhiUnknown = g.fhPhiUnknown; + fhEtaUnknown = g.fhEtaUnknown; + + /* + fhMCElePt = g.fhMCElePt; + fhMCElePhi = g.fhMCElePhi; + fhMCEleEta = g.fhMCEleEta; + */ + + return *this; + +} + +//____________________________________________________________________________ +AliAnaElectron::~AliAnaElectron() +{ + //dtor + +} + + +//________________________________________________________________________ +TList * AliAnaElectron::GetCreateOutputObjects() +{ + // Create histograms to be saved in output file and + // store them in outputContainer + TList * outputContainer = new TList() ; + outputContainer->SetName("ElectronHistos") ; + + Int_t nptbins = GetHistoNPtBins(); + Int_t nphibins = GetHistoNPhiBins(); + Int_t netabins = GetHistoNEtaBins(); + Float_t ptmax = GetHistoPtMax(); + Float_t phimax = GetHistoPhiMax(); + Float_t etamax = GetHistoEtaMax(); + Float_t ptmin = GetHistoPtMin(); + Float_t phimin = GetHistoPhiMin(); + Float_t etamin = GetHistoEtaMin(); + + fh1pOverE = new TH1F("h1pOverE","EMCAL-TRACK matches p/E",100,0.,10.); + fh1dR = new TH1F("h1dR","EMCAL-TRACK matches dR",300, 0.,TMath::Pi()); + fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",200,0.,50.,200,0.,400.); + fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",200,0.,50.,200,0.,400.); + fh2dEtadPhi = new TH2F("h2dEtadPhi","#Delta#eta vs. #Delta#phi for all track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi()); + fh2dEtadPhiMatched = new TH2F("h2dEtadPhiMatched","#Delta#eta vs. #Delta#phi for matched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi()); + fh2dEtadPhiUnmatched = new TH2F("h2dEtadPhiUnmatched","#Delta#eta vs. #Delta#phi for unmatched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi()); + fh2OuterPtVsExtrapPt = new TH2F("h2OuterPtVsExtrapPt","h2OuterPtVsExtrapPt",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); + fh2OuterPhiVsExtrapPhi = new TH2F("h2OuterPhiVsExtrapPhi","h2OuterPhiVsExtrapPhi",nphibins,phimin,phimax,nphibins,phimin,phimax); + fh2OuterEtaVsExtrapEta = new TH2F("h2OuterEtaVsExtrapEta","h2OuterEtaVsExtrapEta",netabins,etamin,etamax,netabins,etamin,etamax); + + fh2TrackPVsClusterE = new TH2F("h2TrackPVsClusterE","h2TrackPVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); + fh2TrackPtVsClusterE = new TH2F("h2TrackPtVsClusterE","h2TrackPtVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); + fh2TrackPhiVsClusterPhi = new TH2F("h2TrackPhiVsClusterPhi","h2TrackPhiVsClusterPhi",nphibins,phimin,phimax,nphibins,phimin,phimax); + fh2TrackEtaVsClusterEta = new TH2F("h2TrackEtaVsClusterEta","h2TrackEtaVsClusterEta",netabins,etamin,etamax,netabins,etamin,etamax); + + outputContainer->Add(fh1pOverE) ; + outputContainer->Add(fh1dR) ; + outputContainer->Add(fh2EledEdx) ; + outputContainer->Add(fh2MatchdEdx) ; + outputContainer->Add(fh2dEtadPhi) ; + outputContainer->Add(fh2dEtadPhiMatched) ; + outputContainer->Add(fh2dEtadPhiUnmatched) ; + outputContainer->Add(fh2OuterPtVsExtrapPt) ; + outputContainer->Add(fh2OuterPhiVsExtrapPhi) ; + outputContainer->Add(fh2OuterEtaVsExtrapEta) ; + outputContainer->Add(fh2TrackPVsClusterE) ; + outputContainer->Add(fh2TrackPtVsClusterE) ; + outputContainer->Add(fh2TrackPhiVsClusterPhi) ; + outputContainer->Add(fh2TrackEtaVsClusterEta) ; + + fhPtElectron = new TH1F("hPtElectron","Electron pT",nptbins,ptmin,ptmax); + fhPhiElectron = new TH2F("hPhiElectron","Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaElectron = new TH2F("hEtaElectron","Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + + outputContainer->Add(fhPtElectron) ; + outputContainer->Add(fhPhiElectron) ; + outputContainer->Add(fhEtaElectron) ; + + if(IsDataMC()){ + + fhPtConversion = new TH1F("hPtConversion","Conversion electron pT",nptbins,ptmin,ptmax); + fhPhiConversion = new TH2F("hPhiConversion","Conversion Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtBottom = new TH1F("hPtBottom","Bottom electron pT",nptbins,ptmin,ptmax); + fhPhiBottom = new TH2F("hPhiBottom","Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaBottom = new TH2F("hEtaBottom","Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtCharm = new TH1F("hPtCharm","Charm electron pT",nptbins,ptmin,ptmax); + fhPhiCharm = new TH2F("hPhiCharm","Charm Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaCharm = new TH2F("hEtaCharm","Charm Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtCFromB = new TH1F("hPtCFromB","Charm from Bottom electron pT",nptbins,ptmin,ptmax); + fhPhiCFromB = new TH2F("hPhiCFromB","Charm from Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaCFromB = new TH2F("hEtaCFromB","Charm from Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtDalitz = new TH1F("hPtDalitz","Dalitz electron pT",nptbins,ptmin,ptmax); + fhPhiDalitz = new TH2F("hPhiDalitz","Dalitz Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaDalitz = new TH2F("hEtaDalitz","Dalitz Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtWDecay = new TH1F("hPtWDecay","W-boson Electron pT",nptbins,ptmin,ptmax); + fhPhiWDecay = new TH2F("hPhiWDecay","W-boson electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaWDecay = new TH2F("hEtaWDecay","W-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtZDecay = new TH1F("hPtZDecay","Z-boson electron pT",nptbins,ptmin,ptmax); + fhPhiZDecay = new TH2F("hPhiZDecay","Z-boson Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaZDecay = new TH2F("hEtaZDecay","Z-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtPrompt = new TH1F("hPtPrompt","Prompt electron pT",nptbins,ptmin,ptmax); + fhPhiPrompt = new TH2F("hPhiPrompt","Prompt Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaPrompt = new TH2F("hEtaPrompt","Prompt Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhPtUnknown = new TH1F("hPtUnknown","Unknown electron pT",nptbins,ptmin,ptmax); + fhPhiUnknown = new TH2F("hPhiUnknown","Unknown Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhEtaUnknown = new TH2F("hEtaUnknown","Unknown Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + + outputContainer->Add(fhPtConversion); + outputContainer->Add(fhPhiConversion); + outputContainer->Add(fhEtaConversion); + outputContainer->Add(fhPtBottom); + outputContainer->Add(fhPhiBottom); + outputContainer->Add(fhEtaBottom); + outputContainer->Add(fhPtCharm); + outputContainer->Add(fhPhiCharm); + outputContainer->Add(fhEtaCharm); + outputContainer->Add(fhPtCFromB); + outputContainer->Add(fhPhiCFromB); + outputContainer->Add(fhEtaCFromB); + outputContainer->Add(fhPtDalitz); + outputContainer->Add(fhPhiDalitz); + outputContainer->Add(fhEtaDalitz); + outputContainer->Add(fhPtWDecay); + outputContainer->Add(fhPhiWDecay); + outputContainer->Add(fhEtaWDecay); + outputContainer->Add(fhPtZDecay); + outputContainer->Add(fhPhiZDecay); + outputContainer->Add(fhEtaZDecay); + outputContainer->Add(fhPtPrompt); + outputContainer->Add(fhPhiPrompt); + outputContainer->Add(fhEtaPrompt); + outputContainer->Add(fhPtUnknown); + outputContainer->Add(fhPhiUnknown); + outputContainer->Add(fhEtaUnknown); + + /* + fhMCElePt = new TH1F("hMCElePt","MC Electron pT",nptbins,ptmin,ptmax); + fhMCElePhi = new TH2F("hMCElePhi","MC Electron phi",nptbins,ptmin,ptmax,nphibins,phimin,phimax); + fhMCEleEta = new TH2F("hMCEleEta","MC Electron eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + + outputContainer->Add(fhMCElePt) ; + outputContainer->Add(fhMCElePhi) ; + outputContainer->Add(fhMCEleEta) ; + */ + + }//Histos with MC + + //Save parameters used for analysis + TString parList ; //this will be list of parameters used for this analysis. + char onePar[255] ; + + sprintf(onePar,"--- AliAnaElectron ---\n") ; + parList+=onePar ; + sprintf(onePar,"fCalorimeter: %s\n",fCalorimeter.Data()) ; + parList+=onePar ; + sprintf(onePar,"fpOverEmin: %f\n",fpOverEmin) ; + parList+=onePar ; + sprintf(onePar,"fpOverEmax: %f\n",fpOverEmax) ; + parList+=onePar ; + sprintf(onePar,"fResidualCut: %f\n",fResidualCut) ; + parList+=onePar ; + + //Get parameters set in base class. + parList += GetBaseParametersList() ; + + //Get parameters set in FidutialCut class (not available yet) + //parlist += GetFidCut()->GetFidCutParametersList() + + TObjString *oString= new TObjString(parList) ; + outputContainer->Add(oString); + + return outputContainer ; + +} + +//____________________________________________________________________________ +void AliAnaElectron::Init() +{ + + //Init + //Do some checks +// if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){ +// printf("AliAnaElectron::Init() - !!ABORT: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n"); +// abort(); +// } +// else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){ +// printf("AliAnaElectron::Init() - !!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n"); +// abort(); +// } + +} + + +//____________________________________________________________________________ +void AliAnaElectron::InitParameters() +{ + + //Initialize the parameters of the analysis. + SetOutputAODClassName("AliAODPWG4Particle"); + SetOutputAODName("PWG4Particle"); + + AddToHistogramsName("AnaElectron_"); + + fCalorimeter = "EMCAL" ; + fpOverEmin = 0.5; + fpOverEmax = 1.5; + fResidualCut = 0.02; + +} + +//__________________________________________________________________ +void AliAnaElectron::MakeAnalysisFillAOD() +{ + // + // Do analysis and fill aods with electron candidates + // These AODs will be used to do subsequent histogram filling + // + // Also fill some QA histograms + // + + //Search for electrons in fCalorimeter + //TRefArray * caloData = new TRefArray(); + //TRefArray * ctsData = new TRefArray(); + cout<<"Event type "<GetInputEvent()->GetName()<GetInputEvent()->GetName(),"AliESDEvent"))) { + printf("AliAnaElectron::MakeAnalysisFillAOD() - !!ABORT: Analysis working only with ESDs!!\n"); + abort(); + } + + //Get vertex for cluster momentum calculation + Double_t vertex[]={0,0,0} ; //vertex ; + if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex); + //Get bfield for track extrapolation to Calorimeter + Double_t bfield = GetReader()->GetInputEvent()->GetMagneticField(); //from V0 finder + Double_t radius = 0.; + + //Get the CTS tracks + //ctsData = GetAODCTS(); + //Select the Calorimeter of the electron + if(fCalorimeter == "PHOS") { + //caloData = GetAODPHOS(); + radius = 425.0; //FIXME + } else if (fCalorimeter == "EMCAL") { + //caloData = GetAODEMCAL(); + radius = 441.0; //[cm] EMCAL radius +13cm FIXME + } + + //if(!(ctsData && caloData) || (ctsData->GetEntriesFast() == 0 || caloData->GetEntriesFast() == 0)) return; + //if(!caloData || caloData->GetEntriesFast() == 0) return; + + //////////////////////////////////////////////// + //Start from tracks and get associated clusters + // + //Note: an alternative method would be to start from clusters and get associated tracks - + //which is better? For electrons, probably tracks-->clusters + //////////////////////////////////////////////// +// for(Int_t itrk = 0; itrk < ctsData->GetEntriesFast(); itrk++){ +// AliAODTrack *track = (AliAODTrack*)ctsData->At(itrk); + + AliESDEvent *esd = (AliESDEvent*) GetReader()->GetInputEvent(); + Int_t nTracks = esd->GetNumberOfTracks() ; + for (Int_t itrk = 0; itrk < nTracks; itrk++) {////////////// track loop + AliESDtrack * track = (AliESDtrack*) esd->GetTrack(itrk) ; + //extrapolate track to Calorimeter + Double_t emcmom[3] = {0.,0.,0.}; + Double_t emcpos[3] = {0.,0.,0.}; + AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam(); + if(!outerparam) continue; + + Bool_t okpos = outerparam->GetXYZAt(radius,bfield,emcpos); + Bool_t okmom = outerparam->GetPxPyPzAt(radius,bfield,emcmom); + if(!(okpos && okmom)) continue; + + TVector3 pos(emcpos[0],emcpos[1],emcpos[2]); + TVector3 mom(emcmom[0],emcmom[1],emcmom[2]); + Double_t tphi = pos.Phi(); + Double_t teta = pos.Eta(); + Double_t tmom = mom.Mag(); + + TLorentzVector mom2(mom,0.); + Bool_t in = GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter) ; + if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, eta %2.2f in fidutial cut %d\n",track->Pt(), track->Phi(), track->Eta(), in); + if(mom.Pt() > GetMinPt() && in) { + + printf("\tExtrapolated pt %2.2f, phi %2.2f, eta %2.2f \n",mom.Pt(),mom.Phi(),mom.Eta()); + fh2OuterPtVsExtrapPt->Fill(mom.Pt(),track->Pt()); + fh2OuterPhiVsExtrapPhi->Fill(mom.Phi(),track->Phi()); + fh2OuterEtaVsExtrapEta->Fill(mom.Eta(),track->Eta()); + + Int_t ntot = esd->GetNumberOfCaloClusters();//caloData->GetEntriesFast(); + Double_t res = 999.; + Double_t pOverE = -999.; + + //For tracks in EMCAL acceptance, pair them with all clusters + //and fill the dEta vs dPhi for these pairs: + for(Int_t iclus = 0; iclus < esd->GetNumberOfCaloClusters(); iclus++) { + AliESDCaloCluster * clus = (AliESDCaloCluster*) esd->GetCaloCluster(iclus); + if(!clus) continue; + if(fCalorimeter == "PHOS" && !clus->IsPHOS()) continue; + if(fCalorimeter == "EMCAL" && !clus->IsEMCAL()) continue; + + Float_t x[3]; + clus->GetPosition(x); + TVector3 cluspos(x[0],x[1],x[2]); + Double_t deta = teta - cluspos.Eta(); + Double_t dphi = tphi - cluspos.Phi(); + if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi(); + if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi(); + if(track->GetEMCALcluster() < -9000) fh2dEtadPhiUnmatched->Fill(deta,dphi); + fh2dEtadPhi->Fill(deta,dphi); + fh2TrackPVsClusterE->Fill(clus->E(),track->P()); + fh2TrackPtVsClusterE->Fill(clus->E(),track->Pt()); + fh2TrackPhiVsClusterPhi->Fill(cluspos.Phi(),mom.Phi()); + fh2TrackEtaVsClusterEta->Fill(cluspos.Eta(),mom.Eta()); + } + + ///////////////////////////////// + //Perform electron cut analysis// + ///////////////////////////////// + Bool_t isElectron = kFALSE; + + Int_t iCluster = track->GetEMCALcluster(); + if(iCluster < -9000) {printf("NOT MATCHED"); continue; }//no match + if(iCluster > ntot) continue; //index out of range; shouldn't happen + if(iCluster < 0 && iCluster > -9000) { //this should only happen in MC events + printf("AliAnaElectron::MakeAnalysisFillAOD() - Track has a fake match: %d\n",iCluster); + continue; + } + AliESDCaloCluster * clus = (AliESDCaloCluster*) esd->GetCaloCluster(iCluster); + if(!clus) continue; + if(fCalorimeter == "PHOS" && !clus->IsPHOS()) continue; + if(fCalorimeter == "EMCAL" && !clus->IsEMCAL()) continue; + + Float_t x[3]; + clus->GetPosition(x); + TVector3 cluspos(x[0],x[1],x[2]); + Double_t deta = teta - cluspos.Eta(); + Double_t dphi = tphi - cluspos.Phi(); + if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi(); + if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi(); + res = sqrt(dphi*dphi + deta*deta); + fh1dR->Fill(res); + fh2dEtadPhiMatched->Fill(deta,dphi); + + if(res < fResidualCut) { + //Good match + fh2MatchdEdx->Fill(track->P(),track->GetTPCsignal()); + + Double_t energy = clus->E(); + if(energy > 0) pOverE = tmom/energy; + fh1pOverE->Fill(pOverE); + + Int_t mult = clus->GetNumberOfDigits(); + // Int_t mcClus = clus->GetLabel(); + AliESDCaloCluster * esdcalo = (AliESDCaloCluster*) esd->GetCaloCluster(clus->GetID()); + Int_t matchIndex = esdcalo->GetTrackMatched(); + + if(matchIndex != itrk) printf("Track and cluster don't agree! track %d, cluster %d",itrk,matchIndex); + if(mult < 2) printf("Single digit cluster."); + + ////////////////////////////// + //Electron cuts happen here!// + ////////////////////////////// + if(pOverE > fpOverEmin && pOverE < fpOverEmax) isElectron = kTRUE; + + } //good matching residual + + /////////////////////////// + //Fill AOD with electrons// + /////////////////////////// + if(isElectron) { + + fh2EledEdx->Fill(track->P(),track->GetTPCsignal()); + + Double_t eMass = 0.511/1000; //mass in GeV + Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass); + AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE); + tr.SetLabel(tr.GetLabel()); + tr.SetCaloLabel(clus->GetID(),-1); //sets the indices of the original caloclusters + tr.SetDetector(fCalorimeter); + tr.SetPdg(11); + + //Play with the MC stack if available + //Check origin of the candidates + if(IsDataMC()){ + + //FIXME: Need to re-think this for track-oriented analysis + tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(clus->GetLabel(),GetMCStack())); + + if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag()); + }//Work with stack also + + AddAODParticle(tr); + + if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg()); + }//electron + }//pt, fiducial selection + }//track loop + + //FIXME: Should we also check from the calocluster side, just in + //case? + + if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs \n"); + +} + +//__________________________________________________________________ +void AliAnaElectron::MakeAnalysisFillHistograms() +{ + //Do analysis and fill histograms + + //Loop on stored AOD electrons + Int_t naod = GetOutputAODBranch()->GetEntriesFast(); + if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod); + + for(Int_t iaod = 0; iaod < naod ; iaod++){ + AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod)); + Int_t pdg = ele->GetPdg(); + + if(GetDebug() > 3) + printf("AliAnaElectron::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ele->GetPdg(),ele->GetTag(), (ele->GetDetector()).Data()) ; + + if(pdg != AliCaloPID::kElectron) continue; + if(ele->GetDetector() != fCalorimeter) continue; + + if(GetDebug() > 1) + printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ele->Pt(),ele->Phi(),ele->Eta()) ; + + //Fill electron histograms + Float_t ptele = ele->Pt(); + Float_t phiele = ele->Phi(); + Float_t etaele = ele->Eta(); + + fhPtElectron ->Fill(ptele); + fhPhiElectron ->Fill(ptele,phiele); + fhEtaElectron ->Fill(ptele,etaele); + + if(IsDataMC()){ + Int_t tag = ele->GetTag(); + if(tag == AliMCAnalysisUtils::kMCConversion){ + fhPtConversion ->Fill(ptele); + fhPhiConversion ->Fill(ptele,phiele); + fhEtaConversion ->Fill(ptele,etaele); + } + else if(tag==AliMCAnalysisUtils::kMCEFromB) + { + fhPtBottom ->Fill(ptele); + fhPhiBottom ->Fill(ptele,phiele); + fhEtaBottom ->Fill(ptele,etaele); + } + else if(tag==AliMCAnalysisUtils::kMCEFromC) + { + fhPtCharm ->Fill(ptele); + fhPhiCharm ->Fill(ptele,phiele); + fhEtaCharm ->Fill(ptele,etaele); + } + else if(tag==AliMCAnalysisUtils::kMCEFromCFromB) + { + fhPtCFromB ->Fill(ptele); + fhPhiCFromB ->Fill(ptele,phiele); + fhEtaCFromB ->Fill(ptele,etaele); + } + else if(tag==AliMCAnalysisUtils::kMCPi0Decay || tag==AliMCAnalysisUtils::kMCEtaDecay || tag==AliMCAnalysisUtils::kMCOtherDecay) + { + fhPtDalitz ->Fill(ptele); + fhPhiDalitz ->Fill(ptele,phiele); + fhEtaDalitz ->Fill(ptele,etaele); + } + else if(tag==AliMCAnalysisUtils::kMCWDecay) + { + fhPtWDecay ->Fill(ptele); + fhPhiWDecay ->Fill(ptele,phiele); + fhEtaWDecay ->Fill(ptele,etaele); + } + else if(tag==AliMCAnalysisUtils::kMCZDecay) + { + fhPtZDecay ->Fill(ptele); + fhPhiZDecay ->Fill(ptele,phiele); + fhEtaZDecay ->Fill(ptele,etaele); + } + else if(tag==AliMCAnalysisUtils::kMCElectron) + { + fhPtPrompt ->Fill(ptele); + fhPhiPrompt ->Fill(ptele,phiele); + fhEtaPrompt ->Fill(ptele,etaele); + } + else{ + fhPtUnknown ->Fill(ptele); + fhPhiUnknown ->Fill(ptele,phiele); + fhEtaUnknown ->Fill(ptele,etaele); + } + }//Histograms with MC + + }// aod loop + + //////////////////////////////////////////////////////// + //Fill histograms of pure MC kinematics from the stack// + //////////////////////////////////////////////////////// + if(IsDataMC()) { + AliStack * stack = GetMCStack() ; + + if(!stack) + printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n"); + + //FIXME: Fill pure kine histograms here + + } + +} + + +//__________________________________________________________________ +void AliAnaElectron::Print(const Option_t * opt) const +{ + //Print some relevant parameters set for the analysis + + if(! opt) + return; + + printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ; + AliAnaPartCorrBaseClass::Print(" "); + + printf("Calorimeter = %s\n", fCalorimeter.Data()) ; + printf("pOverE range = %f - %f\n",fpOverEmin,fpOverEmax); + printf("residual cut = %f\n",fResidualCut); + printf(" \n") ; + +} + +//________________________________________________________________________ +void AliAnaElectron::ReadHistograms(TList* outputList) +{ + // Needed when Terminate is executed in distributed environment + // Refill analysis histograms of this class with corresponding + // histograms in output list. + + // Histograms of this analsys are kept in the same list as other + // analysis, recover the position of + // the first one and then add the next + Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE")); + + //Read histograms, must be in the same order as in + //GetCreateOutputObject. + fh1pOverE = (TH1F *) outputList->At(index); + fh1dR = (TH1F *) outputList->At(index++); + fh2EledEdx = (TH2F *) outputList->At(index++); + fh2MatchdEdx = (TH2F *) outputList->At(index++); + +} + +//__________________________________________________________________ +void AliAnaElectron::Terminate(TList* outputList) +{ + + //Do some plots to end + //Recover histograms from output histograms list, needed for + //distributed analysis. + // ReadHistograms(outputList); + + printf(" AliAnaElectron::Terminate() *** %s Report:", GetName()) ; + +} + diff --git a/PWG4/PartCorrDep/AliAnaElectron.h b/PWG4/PartCorrDep/AliAnaElectron.h new file mode 100755 index 00000000000..393b53387f3 --- /dev/null +++ b/PWG4/PartCorrDep/AliAnaElectron.h @@ -0,0 +1,177 @@ +#ifndef ALIANAELECTRON_H +#define ALIANAELECTRON_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +/* $Id: $ */ + +//_________________________________________________________________________ +// +// Class for the electron identification. +// Clusters from EMCAL matched to tracks are selected +// and kept in the AOD. Few histograms produced. +// + +//-- Author: J.L. Klay (Cal Poly) + +// --- ROOT system --- +class TH2F ; +class TString ; + +// --- ANALYSIS system --- +#include "AliAnaPartCorrBaseClass.h" + +class TList ; + +class AliAnaElectron : public AliAnaPartCorrBaseClass { + +public: + + AliAnaElectron() ; // default ctor + AliAnaElectron(const AliAnaElectron & g) ; // cpy ctor + AliAnaElectron & operator = (const AliAnaElectron & g) ;//cpy assignment + virtual ~AliAnaElectron() ; //virtual dtor + + TList * GetCreateOutputObjects(); + + void Init(); + + void MakeAnalysisFillAOD() ; + + void MakeAnalysisFillHistograms() ; + + void Print(const Option_t * opt)const; + + TString GetCalorimeter() const {return fCalorimeter ; } + Double_t GetpOverEmin() const {return fpOverEmin ; } + Double_t GetpOverEmax() const {return fpOverEmax ; } + + void SetCalorimeter(TString det) {fCalorimeter = det ; } + void SetpOverEmin(Double_t min) {fpOverEmin = min ; } + void SetpOverEmax(Double_t max) {fpOverEmax = max ; } + void SetResidualCut(Double_t cut) {fResidualCut = cut ; } + + void InitParameters(); + + void Terminate(TList * outputList); + void ReadHistograms(TList * outputList); //Fill histograms with + //histograms in ouput list, + //needed in Terminate. + + private: + TString fCalorimeter; //! Which detector? EMCAL or PHOS + Double_t fpOverEmin; //! Minimum p/E value for Electrons + Double_t fpOverEmax; //! Maximum p/E value for Electrons + Double_t fResidualCut; //! Track-cluster matching distance + + //matching checks + TH1F *fh1pOverE; //! p/E for track-cluster matches + TH1F *fh1dR; //! distance between projected track and cluster + TH2F *fh2EledEdx; //! dE/dx vs. momentum for electron candidates + TH2F *fh2MatchdEdx; //! dE/dx vs. momentum for all matches + TH2F *fh2dEtadPhi; //! DeltaEta vs. DeltaPhi of all track/cluster + //! pairs + TH2F *fh2dEtadPhiMatched; //! DeltaEta vs. DeltaPhi of matched + //! track/cluster pairs + TH2F *fh2dEtadPhiUnmatched; //! DeltaEta vs. DeltaPhi of unmatched track/cluster pairs + TH2F *fh2OuterPtVsExtrapPt; + TH2F *fh2OuterPhiVsExtrapPhi; + TH2F *fh2OuterEtaVsExtrapEta; + + TH2F* fh2TrackPVsClusterE; + TH2F* fh2TrackPtVsClusterE; + TH2F* fh2TrackPhiVsClusterPhi; + TH2F* fh2TrackEtaVsClusterEta; + + //Reconstructed + TH1F * fhPtElectron; //! Number of identified electron vs transverse momentum + TH2F * fhPhiElectron; //! Azimuthal angle of identified electron vs transverse momentum + TH2F * fhEtaElectron; //! Pseudorapidity of identified electron vs tranvserse momentum + + TH1F * fhPtConversion; //! Number of conversion electron vs transverse momentum + TH2F * fhPhiConversion; //! Azimuthal angle of conversion electron vs transverse momentum + TH2F * fhEtaConversion; //! Pseudorapidity of conversion electron vs tranvserse momentum + + TH1F * fhPtBottom; //! Number of bottom electron vs transverse momentum + TH2F * fhPhiBottom; //! Azimuthal angle of bottom electron vs transverse momentum + TH2F * fhEtaBottom; //! Pseudorapidity of bottom electron vs tranvserse momentum + + TH1F * fhPtCharm; //! Number of charm electron vs transverse momentum + TH2F * fhPhiCharm; //! Azimuthal angle of charm electron vs transverse momentum + TH2F * fhEtaCharm; //! Pseudorapidity of charm electron vs tranvserse momentum + + TH1F * fhPtCFromB; //! Number of charm from bottom electron vs transverse momentum + TH2F * fhPhiCFromB; //! Azimuthal angle of charm from bottom electron vs transverse momentum + TH2F * fhEtaCFromB; //! Pseudorapidity of charm from bottom electron vs tranvserse momentum + + TH1F * fhPtDalitz; //! Number of dalitz electron vs transverse momentum + TH2F * fhPhiDalitz; //! Azimuthal angle of dalitz electron vs transverse momentum + TH2F * fhEtaDalitz; //! Pseudorapidity of dalitz electron vs tranvserse momentum + + TH1F * fhPtWDecay; //! Number of W-boson electron vs transverse momentum + TH2F * fhPhiWDecay; //! Azimuthal angle of W-boson electron vs transverse momentum + TH2F * fhEtaWDecay; //! Pseudorapidity of W-boson electron vs tranvserse momentum + + TH1F * fhPtZDecay; //! Number of Z-boson electron vs transverse momentum + TH2F * fhPhiZDecay; //! Azimuthal angle of Z-boson electron vs transverse momentum + TH2F * fhEtaZDecay; //! Pseudorapidity of Z-boson electron vs tranvserse momentum + + TH1F * fhPtPrompt; //! Number of prompt electron vs transverse momentum + TH2F * fhPhiPrompt; //! Azimuthal angle of prompt electron vs transverse momentum + TH2F * fhEtaPrompt; //! Pseudorapidity of prompt electron vs tranvserse momentum + + TH1F * fhPtUnknown; //! Number of unknown electron vs transverse momentum + TH2F * fhPhiUnknown; //! Azimuthal angle of unknown electron vs transverse momentum + TH2F * fhEtaUnknown; //! Pseudorapidity of unknown electron vs tranvserse momentum + + /* + //MC + TH1F * fhMCPtElectron; //! pT of MC electrons + TH2F * fhMCPhiElectron; //! Phi of MC electrons + TH2F * fhMCEtaElectron; //! eta of MC electrons + + TH1F * fhMCPtConversion; //! Number of TRACKABLE MC conversion electron vs transverse momentum + TH2F * fhMCPhiConversion; //! Azimuthal angle of TRACKABLE MC conversion electron vs transverse momentum + TH2F * fhMCEtaConversion; //! Pseudorapidity of TRACKABLE MC conversion electron vs tranvserse momentum + + TH1F * fhMCPtBottom; //! Number of MC bottom electron vs transverse momentum + TH2F * fhMCPhiBottom; //! Azimuthal angle of MC bottom electron vs transverse momentum + TH2F * fhMCEtaBottom; //! Pseudorapidity of MC bottom electron vs tranvserse momentum + + TH1F * fhMCPtCharm; //! Number of MC charm electron vs transverse momentum + TH2F * fhMCPhiCharm; //! Azimuthal angle of MC charm electron vs transverse momentum + TH2F * fhMCEtaCharm; //! Pseudorapidity of MC charm electron vs tranvserse momentum + + TH1F * fhMCPtCFromB; //! Number of MC charm from bottom electron vs transverse momentum + TH2F * fhMCPhiCFromB; //! Azimuthal angle of MC charm from bottom electron vs transverse momentum + TH2F * fhMCEtaCFromB; //! Pseudorapidity of MC charm from bottom electron vs tranvserse momentum + + TH1F * fhMCPtDalitz; //! Number of MC dalitz electron vs transverse momentum + TH2F * fhMCPhiDalitz; //! Azimuthal angle of MC dalitz electron vs transverse momentum + TH2F * fhMCEtaDalitz; //! Pseudorapidity of MC dalitz electron vs tranvserse momentum + + TH1F * fhMCPtWDecay; //! Number of MC W-boson electron vs transverse momentum + TH2F * fhMCPhiWDecay; //! Azimuthal angle of MC W-boson electron vs transverse momentum + TH2F * fhMCEtaWDecay; //! Pseudorapidity of MC W-boson electron vs tranvserse momentum + + TH1F * fhMCPtZDecay; //! Number of MC Z-boson electron vs transverse momentum + TH2F * fhMCPhiZDecay; //! Azimuthal angle of MC Z-boson electron vs transverse momentum + TH2F * fhMCEtaZDecay; //! Pseudorapidity of MC Z-boson electron vs tranvserse momentum + + TH1F * fhMCPtPrompt; //! Number of prompt MC electron vs transverse momentum + TH2F * fhMCPhiPrompt; //! Azimuthal angle of prompt MC electron vs transverse momentum + TH2F * fhMCEtaPrompt; //! Pseudorapidity of prompt MC electron vs tranvserse momentum + + TH1F * fhMCPtUnknown; //! Number of unknown MC electron vs transverse momentum + TH2F * fhMCPhiUnknown; //! Azimuthal angle of unknown MC electron vs transverse momentum + TH2F * fhMCEtaUnknown; //! Pseudorapidity of unknown MC electron vs tranvserse momentum + */ + + ClassDef(AliAnaElectron,1) + +} ; + + +#endif//ALIANAELECTRON_H + + + diff --git a/PWG4/libPWG4PartCorrDep.pkg b/PWG4/libPWG4PartCorrDep.pkg index 97ca656939c..615b984982b 100755 --- a/PWG4/libPWG4PartCorrDep.pkg +++ b/PWG4/libPWG4PartCorrDep.pkg @@ -9,7 +9,8 @@ SRCS = PartCorrDep/AliAnaCaloTrigger.cxx \ PartCorrDep/AliAnaParticleJetLeadingConeCorrelation.cxx \ PartCorrDep/AliAnaPhoton.cxx PartCorrDep/AliAnaPi0.cxx \ PartCorrDep/AliAnaPi0EbE.cxx PartCorrDep/AliAnaChargedParticles.cxx \ - PartCorrDep/AliAnaCalorimeterQA.cxx PartCorrDep/AliAnaNeutralMeson.cxx + PartCorrDep/AliAnaCalorimeterQA.cxx PartCorrDep/AliAnaNeutralMeson.cxx \ + PartCorrDep/AliAnaElectron.cxx HDRS:= $(SRCS:.cxx=.h) diff --git a/PWG4/macros/electrons/ConfigAnalysisElectron.C b/PWG4/macros/electrons/ConfigAnalysisElectron.C new file mode 100644 index 00000000000..832e292a049 --- /dev/null +++ b/PWG4/macros/electrons/ConfigAnalysisElectron.C @@ -0,0 +1,89 @@ +/* $Id: $ */ +/* $Log$ */ + +//------------------------------------ +// Configuration macro example: +// +// Do EMCal electron analysis with ESDs +// +//------------------------------------ + +AliAnaPartCorrMaker* ConfigAnalysis() +{ + // + // Configuration goes here + // + printf("======================== \n"); + printf("ConfigAnalysisElectron() \n"); + printf("======================== \n"); + + + //Detector Fidutial Cuts + //AliFidutialCut * fidCut = new AliFidutialCut(); + //fidCut->DoEMCALFidutialCut(kTRUE) ; + + //fidCut->SetSimpleEMCALFidutialCut(0.7,80.,190.); + + //fidCut->Print(""); + + + //----------------------------------------------------------- + // Reader + //----------------------------------------------------------- + AliCaloTrackESDReader *reader = new AliCaloTrackESDReader(); + reader->SetDebug(10);//10 for lots of messages + + //Switch on or off the detectors information that you want + reader->SwitchOffCTS(); + reader->SwitchOffEMCAL(); + reader->SwitchOffEMCALCells(); + reader->SwitchOffPHOS(); + reader->SwitchOffPHOSCells(); + //Min particle pT + //reader->SetEMCALPtMin(1.); + + //reader->SetFidutialCut(fidCut); + reader->Print(""); + + + //--------------------------------------------------------------------- + // Analysis algorithm + //--------------------------------------------------------------------- + + AliAnaElectron *anaelectron = new AliAnaElectron(); + anaelectron->SetDebug(10); //10 for lots of messages + anaelectron->SetCalorimeter("EMCAL"); + anaelectron->SetpOverEmin(0.8); + anaelectron->SetpOverEmax(1.2); + anaelectron->SetResidualCut(0.05); + anaelectron->SwitchOnDataMC(); + anaelectron->SetMinPt(1.); + //Set Histrograms bins and ranges + anaelectron->SetHistoPtRangeAndNBins(0, 50, 100) ; + anaelectron->SetHistoPhiRangeAndNBins(0, TMath::TwoPi(), 100) ; + anaelectron->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ; anaelectron->Print(""); + + //Detector Fidutial Cuts + AliFidutialCut * fidCut2 = new AliFidutialCut(); + fidCut2->DoEMCALFidutialCut(kTRUE) ; + fidCut2->SetSimpleEMCALFidutialCut(0.7,80.,190.); + fidCut2->Print(""); + + + //--------------------------------------------------------------------- + // Set analysis algorithm and reader + //--------------------------------------------------------------------- + maker = new AliAnaPartCorrMaker(); + maker->SetReader(reader);//pointer to reader + maker->AddAnalysis(anaelectron,0); + maker->SetAnaDebug(10) ; + maker->SwitchOnHistogramsMaker() ; + //maker->SwitchOnAODsMaker() ; + + maker->Print(""); + // + printf("======================== \n"); + printf("END ConfigAnalysisElectron() \n"); + printf("======================== \n"); + return maker ; +} diff --git a/PWG4/macros/electrons/GetProofLogs.C b/PWG4/macros/electrons/GetProofLogs.C new file mode 100644 index 00000000000..dc4d431e92d --- /dev/null +++ b/PWG4/macros/electrons/GetProofLogs.C @@ -0,0 +1,7 @@ +{ + gEnv->SetValue("XSec.GSI.DelegProxy","2"); + logs=TProof::Mgr("kread@localhost")->GetSessionLogs(0); + logs->Display(); + //logs->Grep("segmentation violation"); + logs->Save("*","logs.txt"); +} diff --git a/PWG4/macros/electrons/README b/PWG4/macros/electrons/README new file mode 100644 index 00000000000..1516630ccd0 --- /dev/null +++ b/PWG4/macros/electrons/README @@ -0,0 +1,538 @@ +29 May 2009 + +This is a grand unified example of the ALICE Analysis Framework for +EMCal analysis on Local, Local CAF, PROOF, GRID interactive, GRID +batch, and AliEn PLUGIN. + +Notice: Everything is working today. Really! + + +--------------------------------------------------------------------- + +PAR files + +You can skip this PAR section for your first use of the example if you +just want to see how things work. The PAR files in the example were +built with ROOT v5-23-04 and AliRoot v4-16-Rev-11. See VM20090529.pdf +file with updated build instructions for ROOT v5-23-04 and AliROOT +v4-16-Rev-11. If you have write access to $ALICE_ROOT, you can install +AliAnalysisTaskElectron.cxx and make a new PAR file as follows: + +edit $ALICE_ROOT/PWG4/PWG4PartCorrLinkDef.h to include AliAnalysisTaskElectron +edit $ALICE_ROOT/PWG4/libPWG4PartCorr.pkg to include AliAnalysisTaskElectron +edit $ALICE_ROOT/PWG4/CMake_libPWG4PartCorr.txt to include AliAnalysisTaskElectron +cp /AliAnalysisTaskElectron.cxx $ALICE_ROOT/PWG4/PartCorr +cp /AliAnalysisTaskElectron.h $ALICE_ROOT/PWG4/PartCorr + +cd $ALICE_ROOT +make +make PWG4PartCorr.par +mv PWG4PartCorr.par + + +--------------------------------------------------------------------- + +local usage (your laptop, workstation, PDSF account, etc.): + +For PDSF: +ssh -X @pdsf.nersc.gov +cp ~kread/.chos $HOME/.chos (to invoke SciLinux 4.4) +logout and log back in +source /project/projectdirs/alice/alisetup.csh + +All of the files at /afs/cern.ch/user/k/kread/public/all/v416 need to be +copied to your local workarea (except for README and PDF files). + +If not at PDSF, need a valid local matching ROOT and AliROOT installation. +(See /afs/cern.ch/user/k/kread/public/all/v416/VM20090529.pdf) + +cd + +Change this line in anaElectron.C to point to your local AliESDs.root +data files: +char * kInDir = "/afs/cern.ch/user/k/kread/PWG4/data/" + +root -b +.L anaElectron.C +anaElectron(0,"ConfigAnalysisElectron") +.q + + +--------------------------------------------------------------------- + +CAF local usage (including for files on GRID): + +All of the files at /afs/cern.ch/user/k/kread/public/all/v416 need to be +copied to your local workarea. + +Your valid public and private keys need to be properly located in +/afs/cern.ch/user//.globus. + +ssh -X @lxplus.cern.sh +cd +bash +source /afs/cern.ch/alice/caf/caf-lxplus.sh -alien v4-16-Release +enter GRID password + +Edit ESDlist.txt to point to the location of your AliESDs.root files +(including those on the GRID). + +root -b +.L anaElectron.C +anaElectron(1,"ConfigAnalysisElectron") +.q + + +--------------------------------------------------------------------- + +PROOF usage within CERN: + +All of the files at /afs/cern.ch/user/k/kread/public/all/v416 need to be +copied to your local workarea. + +The lines in anaElectron.C to need to look like this (which is the default): + char* myproofname = "alicecaf"; + //char* myproofname = "kread@localhost"; + + +Your valid public and private keys need to be properly located in +$HOME/.globus. + +ssh -X @lxplus.cern.sh (or other computer at CERN) +cd +bash + +You need a valid matching local ROOT and AliROOT installation with ROOT +v5-23-04 or higher. (See +/afs/cern.ch/user/k/kread/public/all/v416/VM20090529.pdf if necessary.) +For instance: +source /afs/cern.ch/alice/caf/caf-lxplus.sh -alien v4-16-Release +enter GRID password + +Edit this line in anaElectron.C to point to the appropriate CAF +dataset: +char * kmydataset = "/COMMON/COMMON/LHC09a4_run8101X" + +To access files in AliEN, first type +root -b +gEnv->SetValue("XSec.GSI.DelegProxy","2") + +root -b +.L anaElectron.C +anaElectron(2,"ConfigAnalysisElectron") +.q + +For instructions concerning how to stage data files to CAF and how +to inspect what data is already staged, see: +http://aliceinfo.cern.ch/Offline/Activities/Analysis/CAF/index.html + +If you do not use the provided PAR files, then you will need to add +the following lines at the beginning of +AliAnalysisTaskElectron::Terminate in AliAnalysisTaskElectron.cxx, + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (mgr->GetAnalysisType()==AliAnalysisManager::kProofAnalysis) return; + +These lines work for all modes of anaElectron.cxx. So, you can install +a revised version of AliAnalysisTaskElectron.cxx at $ALICE_ROOT/PWG4/PartCorr . + +--------------------------------------------------------------------- + +remote PROOF usage (not at CERN): + +This is a GREAT feature, but you will need a local installation with +ROOT v5-23-04 built with --globus, AliROOT v4-16-Rev-nn, and probably +write-access to $ALICE_ROOT. +(See /afs/cern.ch/user/k/kread/public/all/v416/VM20090529.pdf. Remote +PROOF access is not possible without at least ROOT v5-23-02b.) This +requires some understanding of PAR file manipulation. If you have +already been using the example, first type this and be sure not to do it +in the $ALICE_ROOT area: + +cd +/bin/rm -rf ANALYSIS ANALYSISalice AOD ESD STEERBase PWG4PartCorr + +All of the files at /afs/cern.ch/user/k/kread/public/all/v416 need to +be copied to your local workarea. You should only have the (provided) +PAR file versions (not untarred versions) of the AliRoot +subdirectories present in your workarea. + +Modify these lines in anaElectron.C to look like this, replacing +"kread" with your lxplus.cern.ch username: + //char* myproofname = "alicecaf"; + char* myproofname = "kread@localhost"; + +Your valid public and private keys need to be properly located in +$HOME/.globus. + +ssh -L 1093:alicecaf.cern.ch:1093 -fN @lxplus.cern.ch +bash +alien-token-init + +enter GRID password +source /tmp/gclient_env_$UID +export alien_CLOSE_SE="ALICE::CNAF::SE" +export XrdSecGSISRVNAMES="lxfsrd0506.cern.ch" + +Edit this line in anaElectron.C to point to the appropriate CAF +dataset: +char * kmydataset = "/COMMON/COMMON/LHC09a4_run8101X" + +root -b +.L anaElectron.C +anaElectron(2,"ConfigAnalysisElectron") +.q + +For instructions concerning how to stage data files to CAF and how +to inspect what data is already staged, see: +http://aliceinfo.cern.ch/Offline/Activities/Analysis/CAF/index.html . + +The macro will create PAR files in $HOME/.proof/packages. If you have +trouble you can look at what is there and also uncomment appropriate +lines in anaElectron.C to ClearPackages. + +You can see PROOF slave log files by creating a file GetProofLogs.C with these +lines: +{ + gEnv->SetValue("XSec.GSI.DelegProxy","2"); + logs=TProof::Mgr("kread@localhost")->GetSessionLogs(0); + logs->Display(); + //logs->Grep("segmentation violation"); + logs->Save("*","logs.txt"); +} + +and then relaunching root and typing: + +root -b +.x GetProofLogs.C + + +If you do not use the provided PAR files, then you will need to add +the following lines at the beginning of +AliAnalysisTaskElectron::Terminate in AliAnalysisTaskElectron.cxx, + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (mgr->GetAnalysisType()==AliAnalysisManager::kProofAnalysis) return; + +These lines work for all modes of anaElectron.cxx. So, you can install +a revised version of AliAnalysisTaskElectron.cxx at $ALICE_ROOT/PWG4/PartCorr . + +--------------------------------------------------------------------- + +GRID usage: + +Copy all of the files at /afs/cern.ch/user/k/kread/public/all/v416 to your +local workarea. + +Change all instances of username "kread" and to your username in +anaElectron.jdl and mergeElectron.jdl. Change email address in +mergeElectron.jdl. There is no need to change the jobid number inside +mergeElectron.jdl. + +Need a valid local ROOT installation built with --enable-globus. +Need a valid local alien grid shell installed. +(If running locally and not at PDSF or lxplus, see +/afs/cern.ch/user/k/kread/public/all/v416/VM20090529.pdf.) + +Your valid public and private keys need to be properly located in +$HOME/.globus. + +bash +alien-token-init +enter GRID password +source /tmp/gclient_env_$UID +export alien_CLOSE_SE="ALICE::CNAF::SE" + +cd +aliensh (to get to GRID home directory) + +Make the following directories (if not done already): +cd $HOME +mkdir bin +mkdir work3 +mkdir work3/output +mkdir work3/output/merge + +Need to upload files (including par and jdl files) to GRID. + +Can use: +cp file:yourlocalfile.ext gridfile.ext (being sure to spell out the GRID +file name.) +Do not try: +cp file:yourlocalfile.ext . + +Or, exit from aliensh and type locally: +alien_cp localfile.ext alien://alice/cern.ch/user///work3/gridfile.ext + +Do not try to skip the complete GRID file name like this: +alien_cp localfile.ext alien://alice/cern.ch/user///work3/ + +anaElectron.sh needs to be uploaded to the GRID to $HOME/bin. +anaElectron.sh sets the mode to GRID. + +To launch job: +aliensh +cd $HOME/work3 +submit anaElectron.jdl + +To change input data. +aliensh +cd $HOME/work3 +rm mycollect.xml (need to remove before uploading new version) +Use the desired run number and pattern in the find command, perhaps using the -l switch: +find -x mycollect /alice/sim/PDC_08b/LHC08d8/30009/0* AliESDs.root > mycollect.xml +or +find -x mycollect -l 200 /alice/sim/PDC_08b/LHC08d9/40010/* AliESDs.root > mycollect.xml +exit +alien_cp mycollect.xml alien://alice/cern.ch/user///work3/mycollect.xml + +Then, relaunch. + +To monitor grid job: +aliensh +ps +ps -trace +spy stdout +spy stderr +spy anaElectron.log + +If job fails with status E or EV, then to obtain stdout and other output type: +registerOutput + +Available GRID productions monitored via MonaLisa are listed here: +http://pcalimonitor.cern.ch/job_details.jsp + +The choice of $alien_CLOSE_SE is very important. You can see the status of ALICE grid +storage elements here: + +http://pcalimonitor.cern.ch/stats?page=SE/table + + +--------------------------------------------------------------------- + +Plugin usage: + +Copy all of the files at /afs/cern.ch/user/k/kread/public/all/v416 to your +local workarea. + +Change all instances of username "kread" to your username in +anaElectron.jdl and mergeElectron.jdl. Change email address in +mergeElectron.jdl. There is no need to change the jobid number inside +mergeElectron.jdl. + +Need a valid local ROOT installation built with --enable-globus. +Need a valid local alien grid shell installed. +(If running locally and not at PDSF or lxplus, see +/afs/cern.ch/user/k/kread/public/all/v416/VM20090529.pdf.) + +Your valid public and private keys need to be properly located in +$HOME/.globus. + +bash +alien-token-init +enter GRID password +source /tmp/gclient_env_$UID +export alien_CLOSE_SE="ALICE::CNAF::SE" + +Copy new xml file manually to AliEn (after deleting any existing copy in $HOME/work3) by typing +alien_cp mycollect.xml alien://alice/cern.ch/user///work3/mycollect.xml + + +cd +root -b +.L anaElectron.C +anaElectron(4,"ConfigAnalysisElectron") +.q + +Ignore warning messages from plugin about missing myAnalysis.C and +analysis.root which are not needed here. + +Files should be automatically uploaded to grid and job launched. +aliensh is launched in order to monitor job. Exiting aliensh after +jobs are done will automatically invoke merge and copy merged output +histograms to local workarea. + +Note that above we interactively submit the plugin job using mode=4. +However, do not change the mode to 4 inside the file anaElectron.sh +which is executed on each GRID compute node. That needs to remain +mode=3 (GRID). + +See GRID usage section above for more about copying to/from GRID and +monitoring jobs. +Make directories on the GRID (if not done already) as follows: + +aliensh +cd $HOME +mkdir bin +mkdir work3 +mkdir work3/output +mkdir work3/output/merge +exit + + +--------------------------------------------------------------------- + +GRID via AliEn plugin for simpler macro + +An earlier, simpler, separate stand-alone example of EMCal analysis +with the AliEn plugin is presently available at +/afs/cern.ch/user/k/kread/public/plugin. + +Copy the tar file to your local workarea. +Untar the file. +tar -xvf plugin-example.tar + +Need a valid local ROOT installation built with --enable-globus. +Need a valid local alien grid shell installed. +(If running locally and not at PDSF or lxplus, see +/afs/cern.ch/user/k/kread/public/all/v416/VM20090529.pdf.) + +Your valid public and private keys need to be properly located in +$HOME/.globus. + +bash +alien-token-init +enter GRID password +source /tmp/gclient_env_$UID +export alien_CLOSE_SE="ALICE::CNAF::SE" + +cd +root -b +.x runGrid.C +.q + +Edit CreateAlienHandler.C to select different input data available on GRID. +See GRID section above for more about copying to/from GRID and monitoring jobs. +1 June 2009 + +These are the instructions for the anaElctron helper macro when used +with AliRoot v4-17. All of the v416 files and instructions (located at +/afs/cern.ch/user/k/kread/public/all/v416) are appropriate for use +with AliRoot v4-17 except for the following which are different: + +anaElectron.C, anaElectron.jdl, and all of the *.par files. + +The PAR files in this subdirectory were built using ROOT v5-23-04 and +AliRoot v4-17-03. + +The warnings about input and output containers being automatically created are +normal. + +For remote PROOF, follow the usual instructions but set up two tunnels: +ssh -L 1093:alicecaf.cern.ch:1093 -fN @lxplus.cern.ch +ssh -L 11094:alicecaf.cern.ch:11094 -fN @lxplus.cern.ch + +As usual, for PROOF, if you don't use these PAR files, then you will need +to add these lines at the beginning of AliAnalysisTaskElectron::Terminate in +AliAnalysisTaskElectron.cxx: + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (mgr->GetAnalysisType()==AliAnalysisManager::kProofAnalysis) return; + +To build PWG4PartCorrDep.par: + +copy AliAnalysisTaskElectron.cxx and AliAnalysisTaskElectron.h to $ALICE_ROOT/PWG4/PartCorrDep +edit $ALICE_ROOT/PWG4/PWG4PartCorrDepLinkDef.h to include AliAnalysisTaskElectron +edit $ALICE_ROOT/PWG4/libPWG4PartCorrDep.pkg to include AliAnalysisTaskElectron +edit $ALICE_ROOT/PWG4/CMake_libPWG4PartCorrDep.txt to include AliAnalysisTaskElectron + +For PLUGIN mode, we have to work around a new small bug in the AliEn +plugin for v4-17-03. Please, ignore the fact that an irrelevant xml +file (30010.xml) is generated for PLUGIN mode and uploaded to the +grid. This satisfies some (uneccesary) requirements in the plugin for +this release and will be fixed in the future. +June 2009 + +mylauncher.C is a robust file upload and grid submit helper. Just +edit the first few lines (after, of course, defining alien_CLOSE_SE +and doing alien-token-init). Adjust it as needed and you may +prefer it to the plugin. Even for manual grid submitters, it +can really simplify uploading and make it more robust. + +The environment variable alien_CLOSE_SE needs to be defined locally, +and be defined BEFORE you do alien-token-init. (I will change the +other READMEs to reflect this.) Then, and only then, you will find +that it is listed in /tmp/gclient_env_$UID. It needs to be there in +order to be defined later within aliensh. Within aliensh you can see if +things are OK by typing echo $alien_CLOSE_SE. And, of course, outside +of aliensh, you can type echo $alien_CLOSE_SE. + +My recommendation to enforce this is either to put +export alien_CLOSE_SE="ALICE::CNAF::SE" + +in your login profile or initial setup script or to make the following local alias: + +alias myaliti='export alien_CLOSE_SE="ALICE::CNAF::SE";alien-token-init' + +For tcsh/csh, that would of course be + +alias myaliti 'setenv alien_CLOSE_SE "ALICE::CNAF::SE";alien-token-init' + +Then, by just typing "myaliti" (which means my alien token init), +everything will be OK concerning alien_CLOSE_SE and you will never forget +to set it. + +Next, people who use bash will continue to type "source +/tmp/glient_env_$uid" after typing myaliti. + +However, they can to this in their login files + +alias myalige='source /tmp/gclient_env_$UID' + +(That stands for my alien globus environment.) + +Then, they just always type + +myaliti +myalige + +Tcsh/csh users never even need to type "bash" (since they can't source +that bash file in the /tmp area) if they do the following. + +1. Prepare a new file in your util subdirectory and call it +$HOME/util/gclient_env.csh. It should be just like +/tmp/gclient_env_$uid, except remove all the equal signs and change +the word export to setenv everywhere. You only need to do that once. +And you need to make two changes in the new file: + +a. Comment out one line and let the myaliti alias do it: +#setenv alien_CLOSE_SE "ALICE::CNAF::SE" + +b. Replace the corresponding line about the volatile GBBOX_ENVFILE with this sed magic: + +setenv GBBOX_ENVFILE `grep ENVFILE /tmp/gclient_env_$uid | sed s/"export GBBOX_ENVFILE="/""/` + + +2. Then, define: + +alias myalige 'source $HOME/util/gclient_env.csh' + +So, these users never need to type "bash", but also just always type + +myaliti +myalige + + +The symmetry of myaliti, myalige is why I didn't put alien_CLOSE_SE +in my login script. By defining myaliti1, myaliti2, ... myaliti5 +with different, favored SE's in advance, then based on varying storage +element health you can type + +myaliti3 (or 4 or 5 or whatever is best, and avoid mis-typing) +myalige + +I like that. + +Next, note that interactive copying to your alien space does not allow +graceful overwriting. You need to delete any existing copy of a +desired file before you copy a new one to the grid with that name. +(mylauncher.C detects this and does it automatically.) + + + +-------------------------------------------------- + +Note that mylauncher.C opens the way to arbitrary file upload, +download, and copy, and multiple simultaneous, coordinated submits. +That could be useful for systematic analyses involving scores or +hundreds of jobs. diff --git a/PWG4/macros/electrons/anaElectron.C b/PWG4/macros/electrons/anaElectron.C new file mode 100644 index 00000000000..9dfc290baa8 --- /dev/null +++ b/PWG4/macros/electrons/anaElectron.C @@ -0,0 +1,535 @@ +/* $Id: $ */ +//-------------------------------------------------- +// Example macro to do analysis with the +// analysis classes in PWG4PartCorr +// Can be executed with Root and AliRoot +// +// Pay attention to the options and definitions +// set in the lines below +// +// Author : Gustavo Conesa Balbastre (INFN-LNF) +// Modifications by K. Read +// +//------------------------------------------------- +enum anaModes {mLocal, mLocalCAF, mPROOF, mGRID, mPLUGIN}; +//mLocal = 0: Analyze locally files in your computer +//mLocalCAF = 1: Analyze locally CAF files +//mPROOF = 2: Analyze CAF files with PROOF +//mGRID = 3: Analyze files on GRID +//mPLUGIN = 4: Analyze files on GRID with AliEn plugin + +//--------------------------------------------------------------------------- +//Settings to read locally several files, only for "mLocal" mode +//The different values are default, they can be set with environmental +//variables: INDIR, PATTERN, NEVENT, respectively +//NOTE: Must be absolute path. Relative paths don't seem to work +char * kInDir = "/Users/jklay/Projects/LHC/alice/work/ppr/PWG4/v417/data"; +//char * kInDir = "/work/jobs/public/data/"; +char * kPattern = ""; // Data are in files kInDir/kPattern+i +Int_t kFile = 4; // Number of files +//--------------------------------------------------------------------------- +//Collection file for grid analysis +char * kXML = "collection.xml"; +//--------------------------------------------------------------------------- +//Data directory for PROOF analysis +char * kmydataset = "/PWG4/mcosenti/LHC08d10_ppElectronB_Jets#esdTree"; +//char * kmydataset = "/COMMON/COMMON/LHC09a4_run8101X"; +//--------------------------------------------------------------------------- +//Scale histograms from file. Change to kTRUE when xsection file exists +//Put name of file containing xsection +//Put number of events per ESD file +//This is an specific case for normalization of Pythia files. +const Bool_t kGetXSectionFromFileAndScale = kFALSE ; +const char * kXSFileName = "pyxsec.root"; +const Int_t kNumberOfEventsPerFile = 100; +//--------------------------------------------------------------------------- + +const Bool_t kMC = kTRUE; //With real data kMC = kFALSE +const TString kInputData = "ESD";//ESD, AOD, MC +TString kTreeName = "esdTree"; + +void anaElectron(Int_t mode=mLocal, TString configName = "ConfigAnalysisElectron") +{ + // Main + + //-------------------------------------------------------------------- + // Load analysis libraries + // Look at the method below, + // change whatever you need for your analysis case + // ------------------------------------------------------------------ + LoadLibraries(mode) ; + + //------------------------------------------------------------------------------------------------- + //Create chain from ESD and from cross sections files, look below for options. + //------------------------------------------------------------------------------------------------- + if(kInputData == "ESD") kTreeName = "esdTree" ; + else if(kInputData == "AOD") kTreeName = "aodTree" ; + else if (kInputData == "MC") kTreeName = "TE" ; + else { + cout<<"Wrong data type "<SetRunMode("submit"); + //Uncomment the following 3 lines to permit auto xml creation + //plugin->SetGridDataDir("/alice/sim/PDC_08b/LHC08d10/"); //dummy + //plugin->SetDataPattern("AliESDs.root"); //dummy + //plugin->AddRunNumber(30010); //dummy + plugin->AddDataFile("mycollect.xml"); + plugin->SetGridWorkingDir("work3"); + plugin->SetAdditionalLibs("anaElectron.C ConfigAnalysisElectron.C ANALYSIS.par ANALYSISalice.par AOD.par ESD.par PWG4PartCorrBase.par PWG4PartCorrDep.par STEERBase.par"); + plugin->SetJDLName("anaElectron.jdl"); + plugin->SetExecutable("anaElectron.sh"); + plugin->SetOutputFiles("histos.root"); + AliAnalysisGrid *alienHandler = plugin; + if (!alienHandler) return; + + // Connect plug-in to the analysis manager + mgr->SetGridHandler(alienHandler); + } + + // MC handler + if(kMC || kInputData == "MC"){ + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file + mgr->SetMCtruthEventHandler(mcHandler); + if( kInputData == "MC") mgr->SetInputEventHandler(NULL); + } + + // AOD output handler + AliAODHandler* aodoutHandler = new AliAODHandler(); + aodoutHandler->SetOutputFileName("aod.root"); + //aodoutHandler->SetCreateNonStandardAOD(); + mgr->SetOutputEventHandler(aodoutHandler); + + //input + if(kInputData == "ESD"){ + // ESD handler + AliESDInputHandler *esdHandler = new AliESDInputHandler(); + mgr->SetInputEventHandler(esdHandler); + } + if(kInputData == "AOD"){ + // AOD handler + AliAODInputHandler *aodHandler = new AliAODInputHandler(); + mgr->SetInputEventHandler(aodHandler); + } + + mgr->SetDebugLevel(10); // For debugging, do not uncomment if you want no messages. + + //------------------------------------------------------------------------- + //Define task, put here any other task that you want to use. + //------------------------------------------------------------------------- + //AliAnaElectron * taskpwg4 = new AliAnaElectron(); + AliAnalysisTaskParticleCorrelation * taskpwg4 = new AliAnalysisTaskParticleCorrelation ("Particle"); + taskpwg4->SetConfigFileName(configName); //Default name is ConfigAnalysisElectron + + mgr->AddTask(taskpwg4); + + // Create containers for input/output + AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer(); + AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("histos", TList::Class(), + AliAnalysisManager::kOutputContainer, "histos.root"); + + mgr->ConnectInput (taskpwg4, 0, cinput1); + mgr->ConnectOutput (taskpwg4, 0, coutput1 ); + mgr->ConnectOutput (taskpwg4, 1, coutput2 ); + + //------------------------ + //Scaling task + //----------------------- + Int_t nfiles = chainxs->GetEntries(); + //cout<<"Get? "< 0){ + //cout<<"Init AnaScale"<Set(xsection/ntrials/kNumberOfEventsPerFile/nfiles) ; + scale->MakeSumw2(kFALSE);//If you want histograms with error bars set to kTRUE + //scale->SetDebugLevel(2); + mgr->AddTask(scale); + + AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("histosscaled", TList::Class(), + AliAnalysisManager::kOutputContainer, "histosscaled.root"); + mgr->ConnectInput (scale, 0, coutput2); + mgr->ConnectOutput (scale, 0, coutput3 ); + } + + //----------------------- + // Run the analysis + //----------------------- + TString smode = ""; + if (mode==mLocal || mode == mLocalCAF) + smode = "local"; + else if (mode==mPROOF) + smode = "proof"; + else if (mode==mGRID) + smode = "local"; + else if (mode==mPLUGIN) + smode = "grid"; + + //mgr->ResetAnalysis(); + mgr->InitAnalysis(); + mgr->PrintStatus(); + if (mode==mPROOF) + mgr->StartAnalysis(smode.Data(),kmydataset,1500,0); + else if (mode==mPLUGIN) + mgr->StartAnalysis(smode.Data()); + else + mgr->StartAnalysis(smode.Data(),chain); + + cout <<" Analysis ended sucessfully "<< endl ; + + } + else cout << "Chain was not produced ! "<>>>>>>>>>> Local mode <<<<<<<<<<<<<< + //---------------------------------------------------------- + if (mode==mLocal || mode == mLocalCAF || mode == mGRID || mode == mPLUGIN) { + + //-------------------------------------- + // Load the needed libraries most of them already loaded by aliroot + //-------------------------------------- + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libXMLIO.so"); + //-------------------------------------------------------- + // If you want to use already compiled libraries + // in the aliroot distribution + //-------------------------------------------------------- + //gSystem->Load("libSTEERBase"); + //gSystem->Load("libESD"); + //gSystem->Load("libAOD"); + //gSystem->Load("libANALYSIS"); + //gSystem->Load("libANALYSISalice"); + //gSystem->Load("libPWG4PartCorrBase"); + //gSystem->Load("libPWG4PartCorrDep"); + + //-------------------------------------------------------- + //If you want to use root and par files from aliroot + //-------------------------------------------------------- + SetupPar("STEERBase"); + SetupPar("ESD"); + SetupPar("AOD"); + SetupPar("ANALYSIS"); + SetupPar("ANALYSISalice"); + SetupPar("PWG4PartCorrBase"); + SetupPar("PWG4PartCorrDep"); + } + + //--------------------------------------------------------- + // <<<<<<<<<< PROOF mode >>>>>>>>>>>> + //--------------------------------------------------------- + else if (mode==mPROOF) { + // + // Connect to proof + // Put appropriate username here + //char* myproofname = "alicecaf"; + char* myproofname = "jklay@localhost"; + + //TProof::Reset("proof://kread@lxb6046.cern.ch"); + //JLK 28-Jun-2009: Only need to do this occasionally? + //TProof::Reset("jklay@localhost",kTRUE); + + gEnv->SetValue("XSec.GSI.DelegProxy","2"); + //TProof::Mgr(myproofname)->ShowROOTVersions(); + //TProof::Mgr(myproofname)->SetROOTVersion("v5-23-04"); + TProof::Open(myproofname); + + // gProof->ClearPackages(); + // gProof->SetLogLevel(5); + // gProof->ClearPackage("STEERBase"); + // gProof->ClearPackage("ESD"); + // gProof->ClearPackage("AOD"); + // gProof->ClearPackage("ANALYSIS"); + // gProof->ClearPackage("ANALYSISalice"); + // gProof->ClearPackage("PWG4PartCorrBase"); + // gProof->ClearPackage("PWG4PartCorrDep"); + // gProof->ShowEnabledPackages(); + + // Enable the STEERBase Package + gProof->UploadPackage("STEERBase.par"); + gProof->EnablePackage("STEERBase"); + // Enable the ESD Package + gProof->UploadPackage("ESD.par"); + gProof->EnablePackage("ESD"); + // Enable the AOD Package + gProof->UploadPackage("AOD.par"); + gProof->EnablePackage("AOD"); + // Enable the Analysis Package + gProof->UploadPackage("ANALYSIS.par"); + gProof->EnablePackage("ANALYSIS"); + // Enable the Analysis Package + gProof->UploadPackage("ANALYSISalice.par"); + gProof->EnablePackage("ANALYSISalice"); + // Enable PartCorr analysis + gProof->UploadPackage("PWG4PartCorrBase.par"); + gProof->EnablePackage("PWG4PartCorrBase"); + // Enable PartCorr analysis + gProof->UploadPackage("PWG4PartCorrDep.par"); + gProof->EnablePackage("PWG4PartCorrDep"); + + gProof->ShowEnabledPackages(); + } + +} + +void SetupPar(char* pararchivename) +{ + //Load par files, create analysis libraries + //For testing, if par file already decompressed and modified + //classes then do not decompress. + + TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; + TString parpar(Form("%s.par", pararchivename)) ; + //if ( gSystem->AccessPathName(parpar.Data()) ) { + //The lines in this if block are forbidden on GRID. + // gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ; + // TString processline(Form(".! make %s", parpar.Data())) ; + // gROOT->ProcessLine(processline.Data()) ; + // gSystem->ChangeDirectory(cdir) ; + // processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ; + // gROOT->ProcessLine(processline.Data()) ; + //} + if ( gSystem->AccessPathName(pararchivename) ) { + TString processline = Form(".! tar xvzf %s",parpar.Data()) ; + gROOT->ProcessLine(processline.Data()); + } + + TString ocwd = gSystem->WorkingDirectory(); + gSystem->ChangeDirectory(pararchivename); + + // check for BUILD.sh and execute + if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) { + printf("*******************************\n"); + printf("*** Building PAR archive ***\n"); + cout<Exec("PROOF-INF/BUILD.sh")) { + Error("runProcess","Cannot Build the PAR Archive! - Abort!"); + return -1; + } + } + // check for SETUP.C and execute + if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) { + printf("*******************************\n"); + printf("*** Setup PAR archive ***\n"); + cout<Macro("PROOF-INF/SETUP.C"); + } + + gSystem->ChangeDirectory(ocwd.Data()); + printf("Current dir: %s\n", ocwd.Data()); +} + + + +void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){ + //Fills chain with data + TString ocwd = gSystem->WorkingDirectory(); + + //----------------------------------------------------------- + //Analysis of CAF data locally + //----------------------------------------------------------- + if(mode == mLocalCAF){ + // Read the input list of files and add them to the chain + TString line; + ifstream in; + in.open("ESDlist.txt"); + while (in.good()) { + in >> line; + if (line.Length() == 0) continue; + // cout << " line = " << line << endl; + chain->Add(line); + } + } + + //--------------------------------------- + //Local files analysis + //--------------------------------------- + else if(mode == mLocal){ + //If you want to add several ESD files sitting in a common directory INDIR + //Specify as environmental variables the directory (INDIR), the number of files + //to analyze (NEVENT) and the pattern name of the directories with files (PATTERN) + + if(gSystem->Getenv("INDIR")) + kInDir = gSystem->Getenv("INDIR") ; + else cout<<"INDIR not set, use default: "<Getenv("PATTERN")) + kPattern = gSystem->Getenv("PATTERN") ; + else cout<<"PATTERN not set, use default: "<Getenv("NEVENT")) + kFile = atoi(gSystem->Getenv("NEVENT")) ; + else cout<<"NEVENT not set, use default: "<cd(kInDir) ) {//check if ESDs directory exist + printf("%s does not exist\n", kInDir) ; + return ; + } + cout<<"INDIR : "<Get(kTreeName) ) { + printf("++++ Adding %s\n", file) ; + chain->AddFile(file); + chainxs->Add(filexs) ; + } + } + else { + printf("---- Skipping %s\n", file) ; + skipped++ ; + } + } + printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ; + } + else { + TString input = "AliESDs.root" ; + cout<<">>>>>> No list added, take a single file <<<<<<<<< "<AddFile(input); + } + + }// local files analysis + + //------------------------------ + //GRID xml files + //----------------------------- + else if(mode == mGRID){ + //Get colection file. It is specified by the environmental + //variable XML + + if(gSystem->Getenv("XML") ) + kXML = gSystem->Getenv("XML"); + else + sprintf(kXML, "collection.xml") ; + + if (!TFile::Open(kXML)) { + printf("No collection file with name -- %s -- was found\n",kXML); + return ; + } + else cout<<"XML file "<Load("libNetx.so") ; + gSystem->Load("libRAliEn.so"); + TGrid::Connect("alien://") ; + + //Feed Grid with collection file + //TGridCollection * collection = (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML)); + TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML); + if (! collection) { + AliError(Form("%s not found", kXML)) ; + return kFALSE ; + } + TGridResult* result = collection->GetGridResult("",0 ,0); + + // Makes the ESD chain + printf("*** Getting the Chain ***\n"); + for (Int_t index = 0; index < result->GetEntries(); index++) { + TString alienURL = result->GetKey(index, "turl") ; + cout << "================== " << alienURL << endl ; + chain->Add(alienURL) ; + alienURL.ReplaceAll("AliESDs.root",kXSFileName); + chainxs->Add(alienURL) ; + } + }// xml analysis + + gSystem->ChangeDirectory(ocwd.Data()); +} + +//________________________________________________ +void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr) +{ + // Read the PYTHIA statistics from the file pyxsec.root created by + // the function WriteXsection(): + // integrated cross section (xsection) and + // the number of Pyevent() calls (ntrials) + // and calculate the weight per one event xsection/ntrials + // The spectrum calculated by a user should be + // multiplied by this weight, something like this: + // TH1F *userSpectrum ... // book and fill the spectrum + // userSpectrum->Scale(weight) + // + // Yuri Kharlov 19 June 2007 + // Gustavo Conesa 15 April 2008 + Double_t xsection = 0; + UInt_t ntrials = 0; + xs = 0; + ntr = 0; + + Int_t nfiles = tree->GetEntries() ; + if (tree && nfiles > 0) { + tree->SetBranchAddress("xsection",&xsection); + tree->SetBranchAddress("ntrials",&ntrials); + for(Int_t i = 0; i < nfiles; i++){ + tree->GetEntry(i); + xs += xsection ; + ntr += ntrials ; + cout << "xsection " <>>> Empty tree !!!! <<<<< "< anaElectron.log 2>&1 <Load("libNetx.so") ; + gSystem->Load("libRAliEn.so"); + + TGrid::Connect("alien://") ; + if (gGrid && gGrid->IsConnected()) { + TString homedir = gGrid->GetHomeDirectory(); // has a trailing slash + TString workdir = homedir + worksubdir; + if (gGrid->Cd(workdir)) { + + // Upload files listed + if (filelist.Length()) { + arr = filelist.Tokenize(" "); + TObjString *os; + TIter next(arr); + while ((os=(TObjString*)next())) { + Info("Grid Upload", "Copying %s to your AliEn work directory", os->GetString().Data()); + if (FileExists(os->GetString())) gGrid->Rm(os->GetString()); + TFile::Cp(Form("file:%s",os->GetString().Data()), Form("alien://%s/%s", workdir.Data(), os->GetString().Data())); + } + delete arr; + } + + // Upload executable if listed + if (execfilename.Length()) { + filename = Form("%sbin/%s", homedir.Data(), execfilename.Data()); + if (FileExists(filename)) gGrid->Rm(filename); + Info("Grid Upload", "Copying executable file %s to your AliEn bin directory", execfilename.Data()); + TFile::Cp(Form("file:%s",execfilename.Data()), Form("alien://%s", filename.Data())); + } + + // Upload and submit JDL if listed + if (jdlfilename.Length()) { + filename = Form("%s/%s", workdir.Data(), jdlfilename.Data()); + if (FileExists(filename)) gGrid->Rm(filename); + Info("Grid Upload", "Copying JDL file %s to your AliEn work directory", jdlfilename.Data()); + TFile::Cp(Form("file:%s",jdlfilename.Data()), Form("alien://%s", filename.Data())); + + TGridResult *res; + TString jobID = ""; + res = gGrid->Command(Form("submit %s", jdlfilename.Data())); + Info("Launcher:", "Submitting %s ", jdlfilename.Data()); + if (res) { + const char *cjobId = res->GetKey(0,"jobId"); + if (!cjobId) { + Error("Launcher:", "Your JDL %s could not be submitted", jdlfilename.Data()); + return; + } + else { + Info("Launcher:", "Your JDL %s was successfully submitted.\n\n\t\t\t THE JOB ID IS: %s\n", + jdlfilename.Data(), cjobId); + } + delete res; + } + } + + // Launch alien shell + gSystem->Exec("aliensh"); + + } + } +} + +Bool_t FileExists(const char *lfn) const +{ +// Returns true if file exists. + if (!gGrid) { + Error("FileExists", "No connection to grid"); + return kFALSE; + } + TGridResult *res = gGrid->Ls(lfn); + if (!res) return kFALSE; + TMap *map = dynamic_cast(res->At(0)); + if (!map) { + delete res; + return kFALSE; + } + TObjString *objs = dynamic_cast(map->GetValue("name")); + if (!objs || !objs->GetString().Length()) { + delete res; + return kFALSE; + } + delete res; + return kTRUE; +} diff --git a/PWG4/macros/electrons/parmaker b/PWG4/macros/electrons/parmaker new file mode 100755 index 00000000000..b493c5be256 --- /dev/null +++ b/PWG4/macros/electrons/parmaker @@ -0,0 +1,104 @@ +#!/bin/bash + +# +# Author: K. Read +# +# This script makes PAR files without requiring write access to $ALICE_ROOT. +# Execute this script in your local area with write access. +# Usage: +# parmaker ANALYSIS +# parmaker ANALYSIS remote +# parmaker PWG4PartCorrDep +# parmaker PWG4PartCorrDep remote +# +# So far only available for these par files: +# ANALYSIS ANALYSISalice AOD ESD PWG4PartCorrBase PWG4PartCorrDep STEERBase +# + +case $2 in + + "remote") + + case $1 in + "ANALYSIS") + parmaker_input_basedir="ANALYSIS" + parmaker_input_dir_subdir="ANALYSIS" + parmaker_output_dir_subdir="ANALYSIS" + ;; + "ANALYSISalice") + parmaker_input_basedir="ANALYSIS" + parmaker_input_dir_subdir="ANALYSIS" + parmaker_output_dir_subdir="ANALYSISalice" + ;; + "AOD") + parmaker_input_basedir="STEER" + parmaker_input_dir_subdir="STEER" + parmaker_output_dir_subdir="AOD" + ;; + "ESD") + parmaker_input_basedir="STEER" + parmaker_input_dir_subdir="STEER" + parmaker_output_dir_subdir="ESD" + ;; + "PWG4PartCorrBase") + parmaker_input_basedir="PWG4" + parmaker_input_dir_subdir="PWG4/PartCorrBase" + parmaker_output_dir_subdir="PWG4PartCorrBase/PartCorrBase" + ;; + "PWG4PartCorrDep") + parmaker_input_basedir="PWG4" + parmaker_input_dir_subdir="PWG4/PartCorrDep" + parmaker_output_dir_subdir="PWG4PartCorrDep/PartCorrDep" + ;; + "STEERBase") + parmaker_input_basedir="STEER" + parmaker_input_dir_subdir="STEER" + parmaker_output_dir_subdir="STEERBase" + ;; + *) + echo "parmaker: I'm sorry Dave, I'm afraid I can't do that." + exit + esac + + + echo "parmaker to use source $ALICE_ROOT/$parmaker_input_dir_subdir" + if [ -e "$ALICE_ROOT/$parmaker_input_dir_subdir" ] + then + echo "parmaker creating $1.par" + mkdir $1 + if [ $parmaker_input_basedir != $parmaker_input_dir_subdir ] + then + mkdir $parmaker_output_dir_subdir + fi + + list=`grep Ali ${ALICE_ROOT}/${parmaker_input_basedir}/lib${1}.pkg | sed -e 's:.cxx::g' -e 's:SRCS::' -e 's:=::' -e 's:+::' -e 's:\\\::'` + for i in $list; do + cp $ALICE_ROOT/$parmaker_input_basedir/$i.cxx $parmaker_output_dir_subdir + cp $ALICE_ROOT/$parmaker_input_basedir/$i.h $parmaker_output_dir_subdir + done + + mkdir $1/PROOF-INF + cp -r $ALICE_ROOT/$parmaker_input_basedir/PROOF-INF.$1/* $1/PROOF-INF + cp $ALICE_ROOT/$parmaker_input_basedir/${1}LinkDef.h $1 + cp $ALICE_ROOT/$parmaker_input_basedir/lib${1}.pkg $1 + cp $ROOTSYS/test/Makefile.arch $1 + cp $ALICE_ROOT/$parmaker_input_basedir/Makefile $1/Makefiletemp + sed -e 's:include \$(ROOTSYS)\/test\/Makefile.arch:include Makefile.arch:' -e "s:PACKAGE = .*:PACKAGE = ${1}:" $1/Makefiletemp > $1/Makefile + + /bin/rm $1/Makefiletemp + tar cfzh $1.par $1 + /bin/rm -rf $1 + fi + ;; + + *) + echo "parmaker to use local source $1" + if [ -e "$1" ] + then + echo "parmaker creating $1.par" + tar cfzh $1.par $1 + else + echo "local subdirectory $1 not found" + fi + +esac -- 2.43.0