2)New class AliMCAnalysisUtils, extracted method CheckOrigin from AliCaloPID, added some other MC generator utilities
3)Event number can be accessed in the analysis classes with method GetEventNumber
4)AliAnaCaloTrigger derives now from AliAnalysisTaskSE, MC additional information done in AliAnaCaloTriggerMC moved here
set(SRCS
PartCorrBase/AliAODPWG4Particle.cxx PartCorrBase/AliAODPWG4ParticleCorrelation.cxx
PartCorrBase/AliNeutralMesonSelection.cxx PartCorrBase/AliFidutialCut.cxx
- PartCorrBase/AliCaloPID.cxx PartCorrBase/AliIsolationCut.cxx
+ PartCorrBase/AliCaloPID.cxx PartCorrBase/AliMCAnalysisUtils.cxx PartCorrBase/AliIsolationCut.cxx
PartCorrBase/AliAnaScale.cxx PartCorrBase/AliAnaPartCorrMaker.cxx
PartCorrBase/AliAnaPartCorrBaseClass.cxx PartCorrBase/AliAnalysisTaskParticleCorrelation.cxx
PartCorrBase/AliCaloTrackReader.cxx PartCorrBase/AliCaloTrackESDReader.cxx
#pragma link C++ class AliNeutralMesonSelection+;
#pragma link C++ class AliFidutialCut+;
#pragma link C++ class AliCaloPID+;
+#pragma link C++ class AliMCAnalysisUtils+;
#pragma link C++ class AliIsolationCut+;
#pragma link C++ class AliAnaPartCorrMaker+;
#pragma link C++ class AliAnaPartCorrBaseClass+;
#include "AliCaloPID.h"
#include "AliFidutialCut.h"
#include "AliIsolationCut.h"
+#include "AliMCAnalysisUtils.h"
#include "AliNeutralMesonSelection.h"
#include "AliLog.h"
#include "AliAODPWG4ParticleCorrelation.h"
fCaloPID = new AliCaloPID();
fFidCut = new AliFidutialCut();
fIC = new AliIsolationCut();
+ fMCUtils = new AliMCAnalysisUtils();
//Initialize parameters
InitParameters();
fOutputAODName(abc.fOutputAODName), fOutputAODClassName(abc.fOutputAODClassName),
fAODCaloClusters(new TClonesArray(*abc.fAODCaloClusters)),
fAODCaloCells(new AliAODCaloCells(*abc.fAODCaloCells)),
- fCaloPID(abc.fCaloPID), fFidCut(abc.fFidCut), fIC(abc.fIC),fNMS(abc.fNMS),
+ fCaloPID(abc.fCaloPID), fFidCut(abc.fFidCut), fIC(abc.fIC),fMCUtils(abc.fMCUtils), fNMS(abc.fNMS),
fHistoNPtBins(abc.fHistoNPtBins), fHistoPtMax(abc.fHistoPtMax), fHistoPtMin(abc.fHistoPtMin),
fHistoNPhiBins(abc.fHistoNPhiBins), fHistoPhiMax(abc.fHistoPhiMax), fHistoPhiMin(abc.fHistoPhiMin),
fHistoNEtaBins(abc.fHistoNEtaBins), fHistoEtaMax(abc.fHistoEtaMax), fHistoEtaMin(abc.fHistoEtaMin)
fCaloPID = abc.fCaloPID;
fFidCut = abc.fFidCut;
fIC = abc.fIC;
+ fMCUtils = abc.fMCUtils;
fNMS = abc.fNMS;
-
+
fInputAODBranch = new TClonesArray(*abc.fInputAODBranch) ;
fInputAODName = abc.fInputAODName;
fOutputAODBranch = new TClonesArray(*abc.fOutputAODBranch) ;
fOutputAODName = abc.fOutputAODName;
fOutputAODClassName = abc.fOutputAODClassName;
-
fHistoNPtBins = abc.fHistoNPtBins; fHistoPtMax = abc.fHistoPtMax; fHistoPtMin = abc.fHistoPtMin;
fHistoNPhiBins = abc.fHistoNPhiBins; fHistoPhiMax = abc.fHistoPhiMax; fHistoPhiMin = abc.fHistoPhiMin;
fHistoNEtaBins = abc.fHistoNEtaBins; fHistoEtaMax = abc.fHistoEtaMax; fHistoEtaMin = abc.fHistoEtaMin;
delete fAODCaloCells ;
}
- if(fReader) delete fReader ;
+ if(fReader) delete fReader ;
if(fCaloPID) delete fCaloPID ;
- if(fFidCut) delete fFidCut ;
- if(fIC) delete fIC ;
- if(fNMS) delete fNMS ;
+ if(fFidCut) delete fFidCut ;
+ if(fIC) delete fIC ;
+ if(fMCUtils) delete fMCUtils ;
+ if(fNMS) delete fNMS ;
}
}
+//__________________________________________________
+Int_t AliAnaPartCorrBaseClass::GetEventNumber() const {
+ //Get current event number
+
+ return fReader->GetEventNumber() ;
+}
+
//__________________________________________________
AliStack * AliAnaPartCorrBaseClass::GetMCStack() const {
//Get stack pointer from reader
class AliCaloPID ;
class AliFidutialCut ;
class AliIsolationCut ;
+class AliMCAnalysisUtils ;
class AliNeutralMesonSelection ;
class AliStack ;
class AliHeader ;
virtual Int_t GetDebug() const { return fDebug ; }
virtual void SetDebug(Int_t d) { fDebug = d ; }
+ virtual Int_t GetEventNumber() const ;
+
virtual AliCaloTrackReader * GetReader() const {return fReader ; }
virtual void SetReader(AliCaloTrackReader * reader) { fReader = reader ; }
+ virtual void Terminate() {;}
+
//analysis AOD branch
virtual TClonesArray * GetCreateOutputAODBranch() ;
- //{return (new TClonesArray("AliAODPWG4Particle",0)) ;}
virtual TString GetInputAODName() const {return fInputAODName ; }
virtual void SetInputAODName(TString name) { fInputAODName = name; }
virtual TString GetOutputAODName() const {return fOutputAODName ; }
virtual AliHeader* GetMCHeader() const ;
virtual AliGenEventHeader* GetMCGenEventHeader() const ;
-
+ //Analysis helpers classes pointers setters and getters
virtual AliCaloPID * GetCaloPID() const {return fCaloPID ;}
virtual void SetCaloPID(AliCaloPID * pid) { fCaloPID = pid ;}
virtual AliIsolationCut * GetIsolationCut() const {return fIC ;}
virtual void SetIsolationCut(AliIsolationCut * fc) { fIC = fc ;}
+ virtual AliMCAnalysisUtils * GetMCAnalysisUtils() const {return fMCUtils ;}
+ virtual void SetMCAnalysisUtils(AliMCAnalysisUtils * mcutils) { fMCUtils = mcutils ;}
+
virtual AliNeutralMesonSelection * GetNeutralMesonSelection() const {return fNMS ;}
virtual void SetNeutralMesonSelection(AliNeutralMesonSelection * nms) { fNMS = nms ;}
void SetPtCutRange(Double_t ptmin, Double_t ptmax)
{ fMaxPt=ptmax; fMinPt=ptmin;}
-
//Histogrammes setters and getters
virtual void SetHistoPtRangeAndNBins(Float_t min, Float_t max, Int_t n) {
fHistoNPtBins = n ;
TClonesArray* fAODCaloClusters ; //! selected PHOS/EMCAL CaloClusters
AliAODCaloCells * fAODCaloCells ; //! selected PHOS/EMCAL CaloCells
+
+ //Analysis helper classes access pointers
AliCaloPID * fCaloPID; // PID calculation
AliFidutialCut * fFidCut; // Acceptance cuts
AliIsolationCut * fIC; // Isolation cut
+ AliMCAnalysisUtils * fMCUtils; // MonteCarlo Analysis utils
AliNeutralMesonSelection * fNMS; // Neutral Meson Selection
//Histograms binning and range
Float_t fHistoEtaMax ; //Maximum value of eta histogram range
Float_t fHistoEtaMin ; //Minimum value of eta histogram range
- ClassDef(AliAnaPartCorrBaseClass,2)
+ ClassDef(AliAnaPartCorrBaseClass,3)
} ;
fAODBranchList->At(iaod)->Clear();
//Tell the reader to fill the data in the 3 detector lists
- fReader->FillInputEvent();
+ fReader->FillInputEvent(iEntry);
//Loop on analysis algorithms
if(fAnaDebug > 0 ) printf("*** Begin analysis *** \n");
return kTRUE ;
}
+
+//________________________________________________________________________
+void AliAnaPartCorrMaker::Terminate()
+{
+ //Execute Terminate of analysis
+ //Do some final plots.
+
+ for(Int_t iana = 0; iana < fAnalysisContainer->GetEntries(); iana++){
+
+ AliAnaPartCorrBaseClass * ana = ((AliAnaPartCorrBaseClass *) fAnalysisContainer->At(iana)) ;
+ ana->Terminate();
+
+ }//Loop on analysis defined
+}
void SwitchOnAODsMaker() { fMakeAOD = kTRUE ; }
void SwitchOffAODsMaker() { fMakeAOD = kFALSE ; }
+ void Terminate();
+
void AddAnalysis(TObject* ana, Int_t n) {
if ( fAnalysisContainer) fAnalysisContainer->AddAt(ana,n);
else { printf("AnalysisContainer not initialized");
void AliAnalysisTaskParticleCorrelation::Terminate(Option_t */*option*/)
{
// Terminate analysis
+ // Do some plots
//
- AliDebug(1,"Do nothing in Terminate");
- //fAna->Terminate();
+
+ fAna->Terminate();
+
}
#include "AliCaloPID.h"
#include "AliAODCaloCluster.h"
#include "AliAODPWG4Particle.h"
-#include "AliStack.h"
-#include "TParticle.h"
ClassImp(AliCaloPID)
fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) ,
fPHOSNeutralWeight(0.), fPHOSWeightFormula(0),
fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
-fDispCut(0.),fTOFCut(0.), fDebug(-1), fMCGenerator("")
+fDispCut(0.),fTOFCut(0.), fDebug(-1)
{
//Ctor
fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula),
fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula),
fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
-fDebug(pid.fDebug),fMCGenerator(pid.fMCGenerator)
+fDebug(pid.fDebug)
{
// cpy ctor
fDispCut = pid.fDispCut;
fTOFCut = pid.fTOFCut;
fDebug = pid.fDebug;
- fMCGenerator = pid.fMCGenerator;
-
+
return *this;
}
}
-//_________________________________________________________________________
-Int_t AliCaloPID::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) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
-
- 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("AliCaloPID::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("AliCaloPID::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 == 321) 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 == 321) return kMCEtaDecay ;
- else return kMCOtherDecay ;
- }//Decay
- }//HERWIG
- else return kMCUnknown;
- }//Status 1 : Pythia generated
- else if(mStatus == 0){
- if(pPdg ==22 || pPdg ==11) 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){
- printf("Origin electron, pT %f\n",mom->Pt());
-
- if(mStatus == 0) return kMCConversion ;
- else return kMCElectron ;
- }
- else return kMCUnknown;
- }//Good label value
- else{
- if(label < 0 ) printf("AliCaloPID::CheckOrigin: *** bad label or no stack ***: label %d \n", label);
- if(label >= stack->GetNtrack()) printf("AliCaloPID::CheckOrigin: *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
- return kMCUnknown;
- }//Bad label
-
- return kMCUnknown;
-
-}
-
//_______________________________________________________________
void AliCaloPID::InitParameters()
{
fDispCut = 1.5;
fTOFCut = 5.e-9;
fDebug = -1;
- fMCGenerator = "PYTHIA";
}
//_______________________________________________________________
printf("TOF cut = %e\n",fTOFCut);
printf("Dispersion cut = %2.2f\n",fDispCut);
printf("Debug level = %d\n",fDebug);
- printf("MC Generator = %s\n",fMCGenerator.Data());
printf(" \n");
ph->SetChargedBit(ntr>0) ; //Temporary cut, should we evaluate distance?
//Set PID pdg
- ph->SetPdg(GetPdg(calo,cluster->PID(),ph->E()));
-
- if(fDebug > 0){
- printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr);
- printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
- ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
- }
+ ph->SetPdg(GetPdg(calo,cluster->PID(),ph->E()));
+
+ if(fDebug > 0){
+ printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr);
+ printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
+ ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
+ }
}
class AliLog ;
class AliAODCaloCluster;
class AliAODPWG4Particle;
-#include "AliStack.h"
class AliCaloPID : public TObject {
kChargedUnknown=321
};
-
- enum mcTypes {kMCPrompt, kMCFragmentation, kMCISR, kMCPi0Decay, kMCEtaDecay, kMCOtherDecay, kMCPi0, kMCEta, kMCElectron, kMCConversion, kMCUnknown};
+ enum TagType {kPi0Decay, kEtaDecay, kOtherDecay, kConversion, kNoTag = -1};
void InitParameters();
- Int_t CheckOrigin(const Int_t label, AliStack * stack) const ;
Int_t GetPdg(const TString calo, const Double_t * pid, const Float_t energy) 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:
Float_t fEMCALPhotonWeight; //Bayesian PID weight for photons in EMCAL
Float_t fTOFCut; //Cut on TOF, used in PID evaluation
Int_t fDebug; //Debug level
- TString fMCGenerator; // MC geneator used to generate data in simulation
ClassDef(AliCaloPID,3)
} ;
}
//____________________________________________________________________________
-void AliCaloTrackMCReader::FillInputEvent()
-{
- //Create list of particles from EMCAL, PHOS and CTS.
+void AliCaloTrackMCReader::FillInputEvent(Int_t iEntry){
+ //Fill the event counter and input lists that are needed, called by the analysis maker.
+
+ fEventNumber = iEntry;
if(fClonesArrayType == kTParticle){
fAODCTS = new TClonesArray("TParticle",0);
else {AliFatal("Wrong clones type");}
- Int_t indexCh = 0 ;
+ Int_t indexCh = 0 ;
Int_t indexEMCAL = 0 ;
- Int_t indexPHOS = 0 ;
+ Int_t indexPHOS = 0 ;
Int_t iParticle = 0 ;
Double_t charge = 0.;
- for (iParticle=0 ; iParticle < GetStack()->GetNtrack() ; iParticle++) {
+ for (iParticle = 0 ; iParticle < GetStack()->GetNtrack() ; iParticle++) {
TParticle * particle = GetStack()->Particle(iParticle);
TLorentzVector momentum;
Float_t p[3];
void GetVertex(Double_t v[3]) const ;
- void FillInputEvent() ;
+ void FillInputEvent(Int_t iEntry) ;
AliVEvent* GetInputEvent() const {return GetMC();}
void SetInputEvent(TObject* esd, TObject* aod, TObject* mc) ;
//____________________________________________________________________________
AliCaloTrackReader::AliCaloTrackReader() :
- TObject(), fDataType(0), fDebug(0), fFidutialCut(0x0),
+ TObject(), fEventNumber(-1), fDataType(0), fDebug(0),
+ fFidutialCut(0x0),
fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0),
fAODCTS(0x0), fAODEMCAL(0x0), fAODPHOS(0x0),
fEMCALCells(0x0), fPHOSCells(0x0),
//____________________________________________________________________________
AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & g) :
- TObject(g), fDataType(g.fDataType), fDebug(g.fDebug),fFidutialCut(g.fFidutialCut),
+ TObject(g), fEventNumber(g.fEventNumber), fDataType(g.fDataType), fDebug(g.fDebug),
+ fFidutialCut(g.fFidutialCut),
fCTSPtMin(g.fCTSPtMin), fEMCALPtMin(g.fEMCALPtMin),fPHOSPtMin(g.fPHOSPtMin),
fAODCTS(new TClonesArray(*g.fAODCTS)),
fAODEMCAL(new TClonesArray(*g.fAODEMCAL)),
if(&source == this) return *this;
- fDataType = source.fDataType ;
- fDebug = source.fDebug ;
-
+ fDataType = source.fDataType ;
+ fDebug = source.fDebug ;
+ fEventNumber = source.fEventNumber ;
fFidutialCut = source.fFidutialCut;
- fCTSPtMin = source.fCTSPtMin ;
+ fCTSPtMin = source.fCTSPtMin ;
fEMCALPtMin = source.fEMCALPtMin ;
- fPHOSPtMin = source.fPHOSPtMin ;
+ fPHOSPtMin = source.fPHOSPtMin ;
- fAODCTS = new TClonesArray(*source.fAODCTS) ;
- fAODEMCAL = new TClonesArray(*source.fAODEMCAL) ;
- fAODPHOS = new TClonesArray(*source.fAODPHOS) ;
+ fAODCTS = new TClonesArray(*source.fAODCTS) ;
+ fAODEMCAL = new TClonesArray(*source.fAODEMCAL) ;
+ fAODPHOS = new TClonesArray(*source.fAODPHOS) ;
fEMCALCells = new TNamed(*source.fEMCALCells) ;
- fPHOSCells = new TNamed(*source.fPHOSCells) ;
+ fPHOSCells = new TNamed(*source.fPHOSCells) ;
fESD = source.fESD;
- fAOD= source.fAOD;
- fMC = source.fMC;
+ fAOD = source.fAOD;
+ fMC = source.fMC;
- fFillCTS = source.fFillCTS;
- fFillEMCAL = source.fFillEMCAL;
- fFillPHOS = source.fFillPHOS;
+ fFillCTS = source.fFillCTS;
+ fFillEMCAL = source.fFillEMCAL;
+ fFillPHOS = source.fFillPHOS;
fFillEMCALCells = source.fFillEMCALCells;
- fFillPHOSCells = source.fFillPHOSCells;
+ fFillPHOSCells = source.fFillPHOSCells;
return *this;
if(fMC)
return fMC->Header();
else{
- printf("header is not available");
+ printf("Header is not available\n");
return 0x0 ;
}
}
if(fMC)
return fMC->GenEventHeader();
else{
- printf("GenEventHeader is not available");
+ printf("GenEventHeader is not available\n");
return 0x0 ;
}
}
}
//___________________________________________________
-void AliCaloTrackReader::FillInputEvent() {
- //Fill the input lists that are needed, called by the analysis maker.
-
+void AliCaloTrackReader::FillInputEvent(Int_t iEntry) {
+ //Fill the event counter and input lists that are needed, called by the analysis maker.
+
+ fEventNumber = iEntry;
if(fFillCTS) FillInputCTS();
if(fFillEMCAL) FillInputEMCAL();
if(fFillPHOS) FillInputPHOS();
virtual Int_t GetDataType() const { return fDataType ; }
virtual void SetDataType(Int_t data ) { fDataType = data ; }
+ virtual Int_t GetEventNumber() const {return fEventNumber ; }
//Minimum pt setters and getters
virtual Float_t GetEMCALPtMin() const { return fEMCALPtMin ; }
void SwitchOnPHOSCells() {fFillPHOSCells = kTRUE ; }
void SwitchOffPHOSCells() {fFillPHOSCells = kFALSE ; }
- virtual void FillInputEvent() ;
+ virtual void FillInputEvent(Int_t iEntry) ;
virtual void FillInputCTS() {;}
virtual void FillInputEMCAL() {;}
virtual void FillInputPHOS() {;}
virtual void SetInputEvent(TObject* /*esd*/, TObject* /*aod*/, TObject* /*mc*/) {;}
protected:
-
+ Int_t fEventNumber; // Event number
Int_t fDataType ; // Select MC:Kinematics, Data:ESD/AOD, MCData:Both
Int_t fDebug; // Debugging level
AliFidutialCut * fFidutialCut; // Acceptance cuts
-
+
Float_t fCTSPtMin; // pT Threshold on charged particles
Float_t fEMCALPtMin; // pT Threshold on emcal clusters
Float_t fPHOSPtMin; // pT Threshold on phos clusters
--- /dev/null
+/**************************************************************************
+ * 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 <TMath.h>
+#include <TString.h>
+#include <TList.h>
+
+//---- ANALYSIS system ----
+#include "AliLog.h"
+#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) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
+ // 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 == 321) 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 == 321) 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("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("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("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("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("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("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");
+
+}
+
+
--- /dev/null
+#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 <TObject.h>
+class TString ;
+class TList ;
+
+//--- AliRoot system ---
+class AliLog ;
+//#include "AliStack.h"
+class AliStack ;
+//#include "AliGenPythiaEventHeader.h"
+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
+
+
+
// Creates an ntuple for 2x2 and NxN triggers
// Each ntuple connects the maximum trigger amplitudes
// and its positions with reconstructed clusters
+// and if MC stack available, with pt of parent.
//
//*-- Yves Schutz (CERN) & Gustavo Conesa Balbastre (INFN-LNF)
//////////////////////////////////////////////////////////////////////////////
-
-#include <TChain.h>
-#include <TFile.h>
+//Root
#include <TNtuple.h>
#include <TVector3.h>
+//Aliroot
#include "AliAnaCaloTrigger.h"
-#include "AliAnalysisManager.h"
-#include "AliESDEvent.h"
+#include "AliStack.h"
#include "AliLog.h"
#include "AliESDCaloCluster.h"
+#include "AliMCEvent.h"
+#include "AliESDEvent.h"
//______________________________________________________________________________
AliAnaCaloTrigger::AliAnaCaloTrigger() :
- fChain(0),
- fESD(0),
fOutputContainer(0),
fCalorimeter("PHOS"),
fNtTrigger22(0),
//______________________________________________________________________________
AliAnaCaloTrigger::AliAnaCaloTrigger(const char *name) :
- AliAnalysisTask(name,"AnaCaloTrigger"),
- fChain(0),
- fESD(0),
+ AliAnalysisTaskSE(name),
fOutputContainer(0),
fCalorimeter("PHOS"),
fNtTrigger22(0),
{
// Constructor.
- // Input slot #0 works with an Ntuple
- DefineInput(0, TChain::Class());
- // Output slot #0 writes into a TH1 container
- DefineOutput(0, TObjArray::Class()) ;
+ // Output slot
+ DefineOutput(1, TList::Class()) ;
}
+
//____________________________________________________________________________
AliAnaCaloTrigger::AliAnaCaloTrigger(const AliAnaCaloTrigger & ct) :
- AliAnalysisTask(ct), fChain(ct.fChain), fESD(ct.fESD),
+ AliAnalysisTaskSE(ct.GetName()),
fOutputContainer(ct.fOutputContainer), fCalorimeter(ct. fCalorimeter),
fNtTrigger22(ct.fNtTrigger22), fNtTriggerNN(ct.fNtTriggerNN)
{
// cpy ctor
- SetName (ct.GetName()) ;
- SetTitle(ct.GetTitle()) ;
-
+
}
//_________________________________________________________________________
{
// assignment operator
- if(&source == this) return *this;
+ //if(&source == this) return *this;
+ this->~AliAnaCaloTrigger();
+ new(this) AliAnaCaloTrigger(source);
- fChain = source.fChain ;
- fESD = source.fESD ;
fOutputContainer = source.fOutputContainer ;
fCalorimeter = source. fCalorimeter ;
fNtTrigger22 = source.fNtTrigger22 ;
AliAnaCaloTrigger::~AliAnaCaloTrigger()
{
// dtor
- //fOutputContainer->Clear() ;
- //delete fOutputContainer ;
-
+ if(fOutputContainer){
+ fOutputContainer->Clear() ;
+ delete fOutputContainer ;
+ }
}
-//______________________________________________________________________________
-void AliAnaCaloTrigger::ConnectInputData(const Option_t*)
-{
- // Initialisation of branch container and histograms
-
- AliInfo(Form("*** Initialization of %s", GetName())) ;
-
- // Get input data
- fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
- if (!fChain) {
- AliError(Form("Input 0 for %s not found\n", GetName()));
- return ;
- }
-
- fESD = new AliESDEvent();
- fESD->ReadFromTree(fChain);
-
-}
-
//________________________________________________________________________
-void AliAnaCaloTrigger::CreateOutputObjects()
+void AliAnaCaloTrigger::UserCreateOutputObjects()
{
// Create the outputs containers
- OpenFile(0) ;
+ OpenFile(1) ;
// create histograms
- fNtTrigger22 = new TNtuple(fCalorimeter+"trigger22", "Trigger data 2x2 patch", "a22:a220:enMax:phEnMax:eta22:phi22:etaMax:phiMax:phEtaMax:phPhiMax");
- fNtTriggerNN = new TNtuple(fCalorimeter+"triggerNN", "Trigger data NxN patch", "aNN:aNN0:enMax:phEnMax:etaNN:phiNN:etaMax:phiMax:phEtaMax:phPhiMax");
+ fNtTrigger22 = new TNtuple(fCalorimeter+"trigger22", "Trigger data 2x2 patch", "a22:a220:ptGen:enMax:phEnMax:eta22:phi22:etaMax:phiMax:phEtaMax:phPhiMax");
+ fNtTriggerNN = new TNtuple(fCalorimeter+"triggerNN", "Trigger data NxN patch", "aNN:aNN0:ptGen:enMax:phEnMax:etaNN:phiNN:etaMax:phiMax:phEtaMax:phPhiMax");
// create output container
- fOutputContainer = new TObjArray(2) ;
+ fOutputContainer = new TList() ;
fOutputContainer->SetName(GetName()) ;
fOutputContainer->AddAt(fNtTrigger22, 0) ;
}
//______________________________________________________________________________
-void AliAnaCaloTrigger::Exec(Option_t *)
+void AliAnaCaloTrigger::UserExec(Option_t *)
{
- // Processing of one event
-
- Long64_t entry = fChain->GetReadEntry() ;
-
- if (!fESD) {
- AliError("fESD is not connected to the input!") ;
- return ;
- }
-
- if ( !((entry-1)%100) )
- AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
-
- // Get trigger information of fCalorimeter
- TArrayF * triggerAmplitudes = 0x0 ;
- TArrayF * triggerPosition = 0x0 ;
- Int_t numberOfCaloClusters = fESD->GetNumberOfCaloClusters() ;
-
- if(fCalorimeter == "PHOS"){
- triggerAmplitudes = fESD->GetPHOSTriggerAmplitudes();
- triggerPosition = fESD->GetPHOSTriggerPosition();
- }
- else if(fCalorimeter == "EMCAL"){
- triggerAmplitudes = fESD->GetEMCALTriggerAmplitudes();
- triggerPosition = fESD->GetEMCALTriggerPosition();
- }
-
- if( triggerAmplitudes && triggerPosition ){
- // trigger amplitudes
- const Float_t ka22 = static_cast<Float_t>(triggerAmplitudes->At(0)) ;
- const Float_t ka22O = static_cast<Float_t>(triggerAmplitudes->At(1)) ;
- const Float_t kaNN = static_cast<Float_t>(triggerAmplitudes->At(2)) ;
- const Float_t kaNNO = static_cast<Float_t>(triggerAmplitudes->At(3)) ;
-
- // trigger position
- const Float_t kx22 = static_cast<Float_t>(triggerPosition->At(0)) ;
- const Float_t ky22 = static_cast<Float_t>(triggerPosition->At(1)) ;
- const Float_t kz22 = static_cast<Float_t>(triggerPosition->At(2)) ;
- const Float_t kxNN = static_cast<Float_t>(triggerPosition->At(3)) ;
- const Float_t kyNN = static_cast<Float_t>(triggerPosition->At(4)) ;
- const Float_t kzNN = static_cast<Float_t>(triggerPosition->At(5)) ;
-
- Float_t enMax = 0. ;
- Float_t phEnMax = 0. ;
- Float_t etaMax = 0.5 ;
- Float_t phiMax = 0. ;
- Float_t phEtaMax = 0.5 ;
- Float_t phPhiMax = 0. ;
-
- TVector3 vpos22(kx22, ky22, kz22) ;
- TVector3 vposNN(kxNN, kyNN, kzNN) ;
- Float_t eta22 = vpos22.Eta() ;
- Float_t phi22 = vpos22.Phi() * TMath::RadToDeg() + 360. ;
- Float_t etaNN = vposNN.Eta() ;
- Float_t phiNN = vposNN.Phi() * TMath::RadToDeg() + 360. ;
-
- Int_t icaloCluster ;
-
- // loop over the Calorimeters Clusters
-
- for(icaloCluster = 0 ; icaloCluster < numberOfCaloClusters ; icaloCluster++) {
-
- AliESDCaloCluster * cluster = fESD->GetCaloCluster(icaloCluster) ;
-
- if (cluster && ( (fCalorimeter == "PHOS" && cluster->IsPHOS()) ||
- (fCalorimeter == "EMCAL" && cluster->IsEMCAL()))) {
-
- Float_t cluEnergy = cluster->E() ;
- Float_t pos[3] ;
- TVector3 vpos ;
-
- cluster->GetPosition( pos ) ;
-
- if ( cluEnergy > enMax) {
- enMax = cluEnergy ;
- vpos.SetXYZ(pos[0], pos[1], pos[2]) ;
- etaMax = vpos.Eta() ;
- phiMax = vpos.Phi() ;
- }
-
- Double_t * pid = cluster->GetPid() ;
-
- if(pid[AliPID::kPhoton] > 0.9) {
- if ( cluEnergy > phEnMax) {
- phEnMax = cluEnergy ;
- vpos.SetXYZ(pos[0], pos[1], pos[2]) ;
- phEtaMax = vpos.Eta() ;
- phPhiMax = vpos.Phi() ;
+ // Processing of one event
+
+ if ( !((Entry()-1)%100) )
+ printf(" Processing event # %lld\n", Entry()) ;
+ AliESDEvent* esd = (AliESDEvent*)InputEvent();
+
+ //Get MC data, if available
+ AliStack* stack = 0x0;
+ if(MCEvent())
+ stack = MCEvent()->Stack();
+
+ // Get trigger information of fCalorimeter
+ TArrayF * triggerAmplitudes = 0x0 ;
+ TArrayF * triggerPosition = 0x0 ;
+ Int_t numberOfCaloClusters = esd->GetNumberOfCaloClusters() ;
+
+ if(fCalorimeter == "PHOS"){
+ triggerAmplitudes = esd->GetPHOSTriggerAmplitudes();
+ triggerPosition = esd->GetPHOSTriggerPosition();
}
- }
- }//if cluster
-
- fNtTrigger22->Fill(ka22, ka22O, enMax, phEnMax, eta22, phi22, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
- fNtTriggerNN->Fill(kaNN, kaNNO, enMax, phEnMax, etaNN, phiNN, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
- }//CaloCluster loop
-
- }//If trigger arrays filled
-
- PostData(0, fOutputContainer);
-
+ else if(fCalorimeter == "EMCAL"){
+ triggerAmplitudes = esd->GetEMCALTriggerAmplitudes();
+ triggerPosition = esd->GetEMCALTriggerPosition();
+ }
+
+ if( triggerAmplitudes && triggerPosition ){
+ // trigger amplitudes
+ const Float_t ka22 = static_cast<Float_t>(triggerAmplitudes->At(0)) ;
+ const Float_t ka22O = static_cast<Float_t>(triggerAmplitudes->At(1)) ;
+ const Float_t kaNN = static_cast<Float_t>(triggerAmplitudes->At(2)) ;
+ const Float_t kaNNO = static_cast<Float_t>(triggerAmplitudes->At(3)) ;
+
+ // trigger position
+ const Float_t kx22 = static_cast<Float_t>(triggerPosition->At(0)) ;
+ const Float_t ky22 = static_cast<Float_t>(triggerPosition->At(1)) ;
+ const Float_t kz22 = static_cast<Float_t>(triggerPosition->At(2)) ;
+ const Float_t kxNN = static_cast<Float_t>(triggerPosition->At(3)) ;
+ const Float_t kyNN = static_cast<Float_t>(triggerPosition->At(4)) ;
+ const Float_t kzNN = static_cast<Float_t>(triggerPosition->At(5)) ;
+
+ //printf("ka22 %f, ka220 %f, kaNN %f, kaNN0 %f\n",ka22,ka22O,kaNN,kaNNO);
+ //printf("kx22 %f, ky22 %f, kz22 %f, kxNN %f, kyNN %f, kzNN %f \n",kx22,ky22,kz22,kxNN,kyNN,kzNN);
+
+ Float_t enMax = 0. ;
+ Float_t phEnMax = 0. ;
+ Float_t etaMax = 0.5 ;
+ Float_t phiMax = 0. ;
+ Float_t phEtaMax = 0.5 ;
+ Float_t phPhiMax = 0. ;
+
+ TVector3 vpos22(kx22, ky22, kz22) ;
+ TVector3 vposNN(kxNN, kyNN, kzNN) ;
+ Float_t eta22 = vpos22.Eta() ;
+ Float_t phi22 = vpos22.Phi() * TMath::RadToDeg() + 360. ;
+ Float_t etaNN = vposNN.Eta() ;
+ Float_t phiNN = vposNN.Phi() * TMath::RadToDeg() + 360. ;
+
+
+ Int_t icaloCluster = 0 ;
+ Int_t labelmax = -1 ;
+ // loop over the Calorimeters Clusters
+
+ for(icaloCluster = 0 ; icaloCluster < numberOfCaloClusters ; icaloCluster++) {
+
+ AliESDCaloCluster * cluster = esd->GetCaloCluster(icaloCluster) ;
+
+ if (cluster && ( (fCalorimeter == "PHOS" && cluster->IsPHOS()) ||
+ (fCalorimeter == "EMCAL" && cluster->IsEMCAL()))) {
+
+ Float_t cluEnergy = cluster->E() ;
+ Float_t pos[3] ;
+ TVector3 vpos ;
+
+ cluster->GetPosition( pos ) ;
+
+ if ( cluEnergy > enMax) {
+ enMax = cluEnergy ;
+ vpos.SetXYZ(pos[0], pos[1], pos[2]) ;
+ etaMax = vpos.Eta() ;
+ phiMax = vpos.Phi() ;
+ labelmax = cluster->GetLabel();
+ }
+
+ Double_t * pid = cluster->GetPid() ;
+
+ if(pid[AliPID::kPhoton] > 0.9) {
+ if ( cluEnergy > phEnMax) {
+ phEnMax = cluEnergy ;
+ vpos.SetXYZ(pos[0], pos[1], pos[2]) ;
+ phEtaMax = vpos.Eta() ;
+ phPhiMax = vpos.Phi() ;
+ }
+ }
+ }//if cluster
+
+ Float_t ptGen = -1;
+ if(stack && labelmax < stack->GetNtrack() && labelmax >= 0 ){
+ TParticle * particle = stack->Particle(labelmax);
+ ptGen = particle->Energy();
+ }
+
+ fNtTrigger22->Fill(ka22, ka22O, ptGen, enMax, phEnMax, eta22, phi22, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
+ fNtTriggerNN->Fill(kaNN, kaNNO, ptGen, enMax, phEnMax, etaNN, phiNN, etaMax, phiMax * TMath::RadToDeg() + 360., phEtaMax, phPhiMax * TMath::RadToDeg() + 360.);
+
+ }//CaloCluster loop
+
+ }//If trigger arrays filled
+
+ PostData(1, fOutputContainer);
+
}
//______________________________________________________________________________
// Creates an ntuple for 2x2 and NxN triggers
// Each ntuple connects the maximum trigger amplitudes
// and its positions with reconstructed clusters
+// and if MC stack available, with pt of parent.
//
//*-- Yves Schutz (CERN) & Gustavo Conesa Balbastre (INFN-LNF)
//////////////////////////////////////////////////////////////////////////////
-#include "AliAnalysisTask.h"
-class TFile ;
+#include "AliAnalysisTaskSE.h"
class TNtuple ;
-class TH1D ;
-class TH1I ;
-class TChain;
-class AliAnalysisManager ;
class AliESDEvent ;
-class AliAnaCaloTrigger : public AliAnalysisTask {
+class AliAnaCaloTrigger : public AliAnalysisTaskSE {
public:
AliAnaCaloTrigger() ;
AliAnaCaloTrigger & operator=(const AliAnaCaloTrigger& source);
virtual ~AliAnaCaloTrigger() ;
- virtual void Exec(Option_t * opt = "") ;
- virtual void ConnectInputData(Option_t *);
- virtual void CreateOutputObjects();
+ virtual void UserExec(Option_t * opt = "") ;
+ virtual void UserCreateOutputObjects();
// virtual void Terminate(Option_t * opt = "") const ;
TString GetCalorimeter() const {return fCalorimeter ; }
void SetCalorimeter(TString calo) {fCalorimeter = calo ; }
private:
- TChain * fChain ; //!pointer to the analyzed TTree or TChain
- AliESDEvent * fESD ; //! Declaration of leave types
-
- TObjArray * fOutputContainer ; //! output data container
-
+
+ TList * fOutputContainer ; //! output data container
TString fCalorimeter ; // "PHOS" or "EMCAL"
// Histograms
TNtuple * fNtTrigger22 ; //Ntuple with 2x2 max dig amplitude and cluster energy, and positions.
TNtuple * fNtTriggerNN ; //Ntuple with NxN max dig amplitude and cluster energy, and positions.
- ClassDef(AliAnaCaloTrigger, 1); // a photon analysis task
+ ClassDef(AliAnaCaloTrigger, 2); // a trigger analysis task
};
#endif // ALIANACALOTRIGGER_H
//#include "Riostream.h"
#include "TClonesArray.h"
#include "TParticle.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TROOT.h"
//---- AliRoot system ----
#include "AliAnaExample.h"
if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
printf("Example: final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
- if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
- printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
+ // if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
+ //printf("Example: final aod cell entries %d\n", GetAODCaloCells()->GetNumberOfCells());
}
}
}//calo cells container exist
}
}
+
+//__________________________________________________________________
+void AliAnaExample::Terminate()
+{
+
+ //Do some plots to end
+
+ printf(" *** %s Report:", GetName()) ;
+ printf(" pt : %5.3f , RMS : %5.3f \n", fhPt->GetMean(), fhPt->GetRMS() ) ;
+
+ TCanvas * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
+ c->Divide(1, 3);
+
+ c->cd(1) ;
+ gPad->SetLogy();
+ fhPt->SetLineColor(2);
+ fhPt->Draw();
+
+ c->cd(2) ;
+ fhPhi->SetLineColor(2);
+ fhPhi->Draw();
+
+ c->cd(3) ;
+ fhEta->SetLineColor(2);
+ fhEta->Draw();
+
+ c->Print("Example.eps");
+
+ char line[1024] ;
+ sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
+ gROOT->ProcessLine(line);
+ sprintf(line, ".!rm -fR *.eps");
+ gROOT->ProcessLine(line);
+
+ printf("!! All the eps files are in %s.tar.gz !!!", GetName());
+
+}
TString GetDetector() const {return fDetector ;}
void SetDetector( TString calo ) {fDetector = calo; }
+ void Terminate();
+
private:
Int_t fPdg ; //identified particle id
//____________________________________________________________________________
AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation() :
AliAnaPartCorrBaseClass(),
- fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
- fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
+ fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),
+ fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
fhDeltaPhiCharged(0), fhDeltaPhiNeutral(0),
fhDeltaEtaCharged(0), fhDeltaEtaNeutral(0),
fhDeltaPhiChargedPt(0), fhDeltaPhiNeutralPt(0),
AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(const AliAnaParticleHadronCorrelation & g) :
AliAnaPartCorrBaseClass(g),
fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
+ fSelectIsolated(g.fSelectIsolated),
fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral),
fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral),
fhDeltaPhiCharged(g.fhDeltaPhiCharged),
fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
-
+ fSelectIsolated = source.fSelectIsolated ;
+
fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ;
fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ;
fhDeltaPhiCharged = source.fhDeltaPhiCharged ;
SetPtCutRange(2,300);
fDeltaPhiMinCut = 1.5 ;
fDeltaPhiMaxCut = 4.5 ;
-
+ fSelectIsolated = kFALSE;
}
//__________________________________________________________________
printf("pT Hadron < %2.2f\n", GetMaxPt()) ;
printf("Phi trigger particle-Hadron < %3.2f\n", fDeltaPhiMaxCut) ;
printf("Phi trigger particle-Hadron > %3.2f\n", fDeltaPhiMinCut) ;
-
-
+ printf("Isolated Trigger? %d\n", fSelectIsolated) ;
}
//____________________________________________________________________________
printf("Particle %d, In Cluster Refs entries %d\n",iaod, (particle->GetRefClusters())->GetEntriesFast());
}
+ if(OnlyIsolated() && !particle->IsIsolated()) continue;
+
//Make correlation with charged hadrons
if((particle->GetRefTracks())->GetEntriesFast() > 0)
MakeChargedCorrelation(particle, (TSeqCollection*) (particle->GetRefTracks()),kTRUE);
void SetDeltaPhiCutRange(Double_t phimin, Double_t phimax)
{fDeltaPhiMaxCut =phimax; fDeltaPhiMinCut =phimin;}
+ Bool_t OnlyIsolated() const {return fSelectIsolated ; }
+ void SelectIsolated(Bool_t select) {fSelectIsolated = select ; }
+
void InitParameters();
void Print(const Option_t * opt) const;
private:
Double_t fDeltaPhiMaxCut ; // Minimum Delta Phi Gamma-Hadron
- Double_t fDeltaPhiMinCut ; // Maximum Delta Phi Gamma-Hadron
-
+ Double_t fDeltaPhiMinCut ; // Maximum Delta Phi Gamma-Hadron
+ Bool_t fSelectIsolated ; // Select only trigger particles isolated
+
//Histograms
TH2F * fhPhiCharged ; //! Phi distribution of selected charged particles
TH2F * fhPhiNeutral ; //! Phi distribution of selected neutral particles
#include "AliIsolationCut.h"\r
#include "AliNeutralMesonSelection.h"\r
#include "AliAODPWG4ParticleCorrelation.h"\r
-#include "AliCaloPID.h"\r
+#include "AliMCAnalysisUtils.h"\r
\r
ClassImp(AliAnaParticleIsolation)\r
\r
if(IsDataMC()){\r
Int_t tag =aod->GetTag();\r
\r
- if(tag == AliCaloPID::kMCPrompt){\r
+ if(tag == AliMCAnalysisUtils::kMCPrompt){\r
fhPtIsoPrompt ->Fill(ptcluster);\r
fhPhiIsoPrompt ->Fill(ptcluster,phicluster);\r
fhEtaIsoPrompt ->Fill(ptcluster,etacluster);\r
}\r
- else if(tag == AliCaloPID::kMCFragmentation)\r
+ else if(tag == AliMCAnalysisUtils::kMCFragmentation)\r
{\r
fhPtIsoFragmentation ->Fill(ptcluster);\r
fhPhiIsoFragmentation ->Fill(ptcluster,phicluster);\r
fhEtaIsoFragmentation ->Fill(ptcluster,etacluster);\r
}\r
- else if(tag == AliCaloPID::kMCPi0Decay)\r
+ else if(tag == AliMCAnalysisUtils::kMCPi0Decay)\r
{\r
fhPtIsoPi0Decay ->Fill(ptcluster);\r
fhPhiIsoPi0Decay ->Fill(ptcluster,phicluster);\r
fhEtaIsoPi0Decay ->Fill(ptcluster,etacluster);\r
}\r
- else if(tag == AliCaloPID::kMCEtaDecay || tag == AliCaloPID::kMCOtherDecay)\r
+ else if(tag == AliMCAnalysisUtils::kMCEtaDecay || tag == AliMCAnalysisUtils::kMCOtherDecay)\r
{\r
fhPtIsoOtherDecay ->Fill(ptcluster);\r
fhPhiIsoOtherDecay ->Fill(ptcluster,phicluster);\r
fhEtaIsoOtherDecay ->Fill(ptcluster,etacluster);\r
}\r
- else if(tag == AliCaloPID::kMCConversion)\r
+ else if(tag == AliMCAnalysisUtils::kMCConversion)\r
{\r
fhPtIsoConversion ->Fill(ptcluster);\r
fhPhiIsoConversion ->Fill(ptcluster,phicluster);\r
if(n[icone][ipt] == 0) {\r
fhPtThresIsolated[icone][ipt]->Fill(ptC);\r
if(IsDataMC()){\r
- if(tag == AliCaloPID::kMCPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCOtherDecay || tag == AliCaloPID::kMCEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;\r
+ if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;\r
else fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC) ;\r
}\r
}\r
if(nfrac[icone][ipt] == 0) {\r
fhPtFracIsolated[icone][ipt]->Fill(ptC);\r
if(IsDataMC()){\r
- if(tag == AliCaloPID::kMCPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
- else if(tag == AliCaloPID::kMCOtherDecay || tag == AliCaloPID::kMCEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;\r
+ if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC) ;\r
else fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC) ;\r
}\r
}\r
//Sum in cone histograms\r
fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;\r
if(IsDataMC()){\r
- if(tag == AliCaloPID::kMCPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;\r
- else if(tag == AliCaloPID::kMCConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;\r
- else if(tag == AliCaloPID::kMCFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;\r
- else if(tag == AliCaloPID::kMCPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;\r
- else if(tag == AliCaloPID::kMCOtherDecay || tag == AliCaloPID::kMCEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;\r
+ if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;\r
+ else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;\r
else fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;\r
}\r
\r
#include "AliAnaPhoton.h"
#include "AliCaloTrackReader.h"
#include "AliCaloPID.h"
+#include "AliMCAnalysisUtils.h"
+#include "AliStack.h"
ClassImp(AliAnaPhoton)
//MC
fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
+ fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
//MC
fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
+ fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR),
fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
- fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
+ fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
+
{
// cpy ctor
fMinDist2 = g.fMinDist2;
fMinDist3 = g.fMinDist3;
- fhPtPhoton = g.fhPtPhoton ;
+ fhPtPhoton = g.fhPtPhoton ;
fhPhiPhoton = g.fhPhiPhoton ;
fhEtaPhoton = g.fhEtaPhoton ;
- fhPtPrompt = g.fhPtPrompt;
+ fhPtPrompt = g.fhPtPrompt;
fhPhiPrompt = g.fhPhiPrompt;
fhEtaPrompt = g.fhEtaPrompt;
- fhPtFragmentation = g.fhPtFragmentation;
+ fhPtFragmentation = g.fhPtFragmentation;
fhPhiFragmentation = g.fhPhiFragmentation;
fhEtaFragmentation = g.fhEtaFragmentation;
- fhPtPi0Decay = g.fhPtPi0Decay;
+ fhPtISR = g.fhPtISR;
+ fhPhiISR = g.fhPhiISR;
+ fhEtaISR = g.fhEtaISR;
+ fhPtPi0Decay = g.fhPtPi0Decay;
fhPhiPi0Decay = g.fhPhiPi0Decay;
fhEtaPi0Decay = g.fhEtaPi0Decay;
- fhPtOtherDecay = g.fhPtOtherDecay;
+ fhPtOtherDecay = g.fhPtOtherDecay;
fhPhiOtherDecay = g.fhPhiOtherDecay;
fhEtaOtherDecay = g.fhEtaOtherDecay;
- fhPtConversion = g. fhPtConversion;
+ fhPtConversion = g. fhPtConversion;
fhPhiConversion = g.fhPhiConversion;
fhEtaConversion = g.fhEtaConversion;
- fhPtUnknown = g.fhPtUnknown;
+ fhPtUnknown = g.fhPtUnknown;
fhPhiUnknown = g.fhPhiUnknown;
fhEtaUnknown = g.fhEtaUnknown;
fhEtaPrompt->SetYTitle("#eta");
fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
outputContainer->Add(fhEtaPrompt) ;
-
+
fhPtFragmentation = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
fhPtFragmentation->SetYTitle("N");
fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
outputContainer->Add(fhEtaFragmentation) ;
+ fhPtISR = new TH1F("hPtISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
+ fhPtISR->SetYTitle("N");
+ fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
+ outputContainer->Add(fhPtISR) ;
+
+ fhPhiISR = new TH2F
+ ("hPhiISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ fhPhiISR->SetYTitle("#phi");
+ fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
+ outputContainer->Add(fhPhiISR) ;
+
+ fhEtaISR = new TH2F
+ ("hEtaISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+ fhEtaISR->SetYTitle("#eta");
+ fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
+ outputContainer->Add(fhEtaISR) ;
+
fhPtPi0Decay = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
fhPtPi0Decay->SetYTitle("N");
fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
fhEtaConversion->SetYTitle("#eta");
fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
outputContainer->Add(fhEtaConversion) ;
-
+
fhPtUnknown = new TH1F("hPtUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
fhPtUnknown->SetYTitle("N");
fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
fhEtaUnknown->SetYTitle("#eta");
fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
outputContainer->Add(fhEtaUnknown) ;
+
}//Histos with MC
//Save parameters used for analysis
}
+//____________________________________________________________________________
+void AliAnaPhoton::Init()
+{
+
+ //Init
+ //Do some checks
+ if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
+ printf("!!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("!!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
+ abort();
+ }
+
+}
+
+
//____________________________________________________________________________
void AliAnaPhoton::InitParameters()
{
//Play with the MC stack if available
//Check origin of the candidates
if(IsDataMC()){
- aodph.SetTag(GetCaloPID()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
+ aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
if(GetDebug() > 0) printf("AliAnaPhoton::FillAOD: Origin of candidate %d\n",aodph.GetTag());
}//Work with stack also
fhPtPhoton ->Fill(ptcluster);
fhPhiPhoton ->Fill(ptcluster,phicluster);
fhEtaPhoton ->Fill(ptcluster,etacluster);
-
+
if(IsDataMC()){
Int_t tag =ph->GetTag();
- if(tag == AliCaloPID::kMCPrompt){
- fhPtPrompt ->Fill(ptcluster);
+ if(tag == AliMCAnalysisUtils::kMCPrompt){
+ fhPtPrompt ->Fill(ptcluster);
fhPhiPrompt ->Fill(ptcluster,phicluster);
fhEtaPrompt ->Fill(ptcluster,etacluster);
}
- else if(tag==AliCaloPID::kMCFragmentation)
+ else if(tag==AliMCAnalysisUtils::kMCFragmentation)
{
- fhPtFragmentation ->Fill(ptcluster);
+ fhPtFragmentation ->Fill(ptcluster);
fhPhiFragmentation ->Fill(ptcluster,phicluster);
fhEtaFragmentation ->Fill(ptcluster,etacluster);
}
- else if(tag==AliCaloPID::kMCPi0Decay)
+ else if(tag==AliMCAnalysisUtils::kMCISR)
{
- fhPtPi0Decay ->Fill(ptcluster);
+ fhPtISR ->Fill(ptcluster);
+ fhPhiISR ->Fill(ptcluster,phicluster);
+ fhEtaISR ->Fill(ptcluster,etacluster);
+ }
+ else if(tag==AliMCAnalysisUtils::kMCPi0Decay)
+ {
+ fhPtPi0Decay ->Fill(ptcluster);
fhPhiPi0Decay ->Fill(ptcluster,phicluster);
fhEtaPi0Decay ->Fill(ptcluster,etacluster);
}
- else if(tag==AliCaloPID::kMCEtaDecay || tag==AliCaloPID::kMCOtherDecay)
+ else if(tag==AliMCAnalysisUtils::kMCEtaDecay || tag==AliMCAnalysisUtils::kMCOtherDecay)
{
- fhPtOtherDecay ->Fill(ptcluster);
+ fhPtOtherDecay ->Fill(ptcluster);
fhPhiOtherDecay ->Fill(ptcluster,phicluster);
fhEtaOtherDecay ->Fill(ptcluster,etacluster);
}
- else if(tag==AliCaloPID::kMCConversion)
+ else if(tag==AliMCAnalysisUtils::kMCConversion)
{
- fhPtConversion ->Fill(ptcluster);
+ fhPtConversion ->Fill(ptcluster);
fhPhiConversion ->Fill(ptcluster,phicluster);
fhEtaConversion ->Fill(ptcluster,etacluster);
+
}
else{
- fhPtUnknown ->Fill(ptcluster);
+ fhPtUnknown ->Fill(ptcluster);
fhPhiUnknown ->Fill(ptcluster,phicluster);
fhEtaUnknown ->Fill(ptcluster,etacluster);
}
TList * GetCreateOutputObjects();
+ void Init();
+
void MakeAnalysisFillAOD() ;
void MakeAnalysisFillHistograms() ;
TH2F * fhPhiFragmentation; //! Phi of identified fragmentation gamma
TH2F * fhEtaFragmentation; //! eta of identified fragmentation gamma
+ TH1F * fhPtISR; //! Number of identified initial state radiation gamma
+ TH2F * fhPhiISR; //! Phi of identified initial state radiation gamma
+ TH2F * fhEtaISR; //! eta of identified initial state radiation gamma
+
TH1F * fhPtPi0Decay; //! Number of identified Pi0Decay gamma
TH2F * fhPhiPi0Decay; //! Phi of identified Pi0Decay gamma
TH2F * fhEtaPi0Decay; //! eta of identified Pi0Decay gamma
TH2F * fhPhiUnknown; //! Phi of identified Unknown gamma
TH2F * fhEtaUnknown; //! eta of identified Unknown gamma
- ClassDef(AliAnaPhoton,1)
+ ClassDef(AliAnaPhoton,2)
+
} ;
// --- ROOT system ---
#include "TH3.h"
//#include "Riostream.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TROOT.h"
//---- AliRoot system ----
#include "AliAnaPi0.h"
#include "AliCaloTrackReader.h"
#include "AliCaloPID.h"
+#include "AliStack.h"
#ifdef __PHOSGEO__
#include "AliPHOSGeoUtils.h"
//____________________________________________________________________________________________________________________________________________________
-void AliAnaPi0:: MakeAnalysisFillHistograms()
+void AliAnaPi0::MakeAnalysisFillHistograms()
{
//Process one event and extract photons from AOD branch
// filled with AliAnaPhoton and fill histos with invariant mass
+//____________________________________________________________________________________________________________________________________________________
+void AliAnaPi0::Terminate()
+{
+ //Do some calculations and plots from the final histograms.
+
+ printf(" *** %s Terminate:", GetName()) ;
+ printf(" Mgg Real : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(), fhRe1[0]->GetRMS() ) ;
+
+ TCanvas * cIM = new TCanvas("cIM", "", 400, 10, 600, 700) ;
+ cIM->Divide(2, 2);
+
+ cIM->cd(1) ;
+ //gPad->SetLogy();
+ TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ();
+ hIMAllPt->SetLineColor(2);
+ hIMAllPt->SetTitle("No cut on p_{T, #gamma#gamma} ");
+ hIMAllPt->Draw();
+
+ cIM->cd(2) ;
+ TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone("IMPt5");
+ hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
+ TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D("z");
+ hIMPt5->SetLineColor(2);
+ hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
+ hIMPt5->Draw();
+
+ cIM->cd(3) ;
+ TH3F * hRe1Pt10 = (TH3F*)fhRe1[0]->Clone("IMPt10");
+ hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
+ TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D("z");
+ hIMPt10->SetLineColor(2);
+ hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
+ hIMPt10->Draw();
+
+ cIM->cd(4) ;
+ TH3F * hRe1Pt20 = (TH3F*)fhRe1[0]->Clone("IMPt20");
+ hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
+ TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D("z");
+ hIMPt20->SetLineColor(2);
+ hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
+ hIMPt20->Draw();
+
+ cIM->Print("Mgg.eps");
+
+
+ TCanvas * cPt = new TCanvas("cPt", "", 400, 10, 600, 700) ;
+ cPt->Divide(2, 2);
+
+ cPt->cd(1) ;
+ //gPad->SetLogy();
+ TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
+ hPt->SetLineColor(2);
+ hPt->SetTitle("No cut on M_{#gamma#gamma} ");
+ hPt->Draw();
+
+ cPt->cd(2) ;
+ TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone("PtIM1");
+ hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
+ TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
+ hPtIM1->SetLineColor(2);
+ hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
+ hPtIM1->Draw();
+
+ cPt->cd(3) ;
+ TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone("PtIM2");
+ hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
+ TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
+ hPtIM2->SetLineColor(2);
+ hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
+ hPtIM2->Draw();
+
+ cPt->cd(4) ;
+ TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone("PtIM3");
+ hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
+ TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
+ hPtIM3->SetLineColor(2);
+ hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
+ hPtIM3->Draw();
+
+ cPt->Print("Pt.eps");
+
+ char line[1024] ;
+ sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
+ gROOT->ProcessLine(line);
+ sprintf(line, ".!rm -fR *.eps");
+ gROOT->ProcessLine(line);
+
+ printf("!! All the eps files are in %s.tar.gz !!!", GetName());
+}
#include "AliAnaPartCorrBaseClass.h"
#ifdef __PHOSGEO__
- class AliPHOSGeoUtils;
+ class AliPHOSGeoUtils;
#endif
class AliAnaPi0 : public AliAnaPartCorrBaseClass {
-
-public:
-
- AliAnaPi0() ; // default ctor
- AliAnaPi0(const char *name) ; // default ctor
- AliAnaPi0(const AliAnaPi0 & g) ; // cpy ctor
- AliAnaPi0 & operator = (const AliAnaPi0 & api0) ;//cpy assignment
- virtual ~AliAnaPi0() ;//virtual dtor
-
- TList * GetCreateOutputObjects();
-
- void Print(const Option_t * opt) const;
-
- void Init();
- void InitParameters();
-
- //void MakeAnalysisFillAOD() {;} //Not needed
- void MakeAnalysisFillHistograms();
-
-// void SetBadRunsList(){;} ; //Set list of runs which can be used for this analysis
- //To be defined in future.
- void SetEtalonHisto(TH3D * h);//Provide etalon of binning for histograms
-
- //Setters for parameters of event buffers
- void SetNCentrBin(Int_t n=5){fNCentrBin=n ;} //number of bins in centrality
- void SetNZvertBin(Int_t n=5){fNZvertBin=n ;} //number of bins for vertex position
- void SetNRPBin(Int_t n=6) {fNrpBin=n ;} //number of bins in reaction plain
- void SetNMaxEvMix(Int_t n=20){fNmaxMixEv=n ;}//Maximal number of events for mixing
-
- //Setters for event selection
- void SetZvertexCut(Float_t zcut=40.){fZvtxCut=zcut ;} //cut on vertex position
-
- TString GetCalorimeter() const {return fCalorimeter ; }
- void SetCalorimeter(TString det) {fCalorimeter = det ; }
-
-private:
- Bool_t IsBadRun(Int_t /*iRun*/) const {return kFALSE;} //Tests if this run bad according to private list
-
-private:
-
- Int_t fNCentrBin ; //number of bins in event container for centrality
- Int_t fNZvertBin ; //number of bins in event container for vertex position
- Int_t fNrpBin ; //number of bins in event container for reaction plain
- Int_t fNPID ; //Number of possible PID combinations
- Int_t fNmaxMixEv ; //Maximal number of events stored in buffer for mixing
- Float_t fZvtxCut ; //Cut on vertex position
- TString fCalorimeter ; //Select Calorimeter for IM
-
- TList ** fEventsList ; //! containers for photons in stored events
-
- //Histograms
-
- TH3D * fhEtalon ; //Etalon histo, all distributions will have same binning as this one
-
- TH3D ** fhRe1 ; //!REAL two-photon invariant mass distribution for different centralities and PID
- TH3D ** fhMi1 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
- TH3D ** fhRe2 ; //!REAL two-photon invariant mass distribution for different centralities and PID
- TH3D ** fhMi2 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
- TH3D ** fhRe3 ; //!REAL two-photon invariant mass distribution for different centralities and PID
- TH3D ** fhMi3 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
- TH3D * fhEvents; //!Number of events per centrality, RP, zbin
-
- //Acceptance
- TH1D * fhPrimPt ; //! Spectrum of Primary
- TH1D * fhPrimAccPt ; //! Spectrum of primary with accepted daughters
- TH1D * fhPrimY ; //! Rapidity distribution of primary particles
- TH1D * fhPrimAccY ; //! Rapidity distribution of primary with accepted daughters
- TH1D * fhPrimPhi ; //! Azimutal distribution of primary particles
- TH1D * fhPrimAccPhi; //! Azimutal distribution of primary with accepted daughters
-
+
+ public:
+
+ AliAnaPi0() ; // default ctor
+ AliAnaPi0(const char *name) ; // default ctor
+ AliAnaPi0(const AliAnaPi0 & g) ; // cpy ctor
+ AliAnaPi0 & operator = (const AliAnaPi0 & api0) ;//cpy assignment
+ virtual ~AliAnaPi0() ;//virtual dtor
+
+ TList * GetCreateOutputObjects();
+
+ void Print(const Option_t * opt) const;
+
+ void Init();
+ void InitParameters();
+
+ //void MakeAnalysisFillAOD() {;} //Not needed
+ void MakeAnalysisFillHistograms();
+
+ // void SetBadRunsList(){;} ; //Set list of runs which can be used for this analysis
+ //To be defined in future.
+ void SetEtalonHisto(TH3D * h);//Provide etalon of binning for histograms
+
+ //Setters for parameters of event buffers
+ void SetNCentrBin(Int_t n=5){fNCentrBin=n ;} //number of bins in centrality
+ void SetNZvertBin(Int_t n=5){fNZvertBin=n ;} //number of bins for vertex position
+ void SetNRPBin(Int_t n=6) {fNrpBin=n ;} //number of bins in reaction plain
+ void SetNMaxEvMix(Int_t n=20){fNmaxMixEv=n ;}//Maximal number of events for mixing
+
+ //Setters for event selection
+ void SetZvertexCut(Float_t zcut=40.){fZvtxCut=zcut ;} //cut on vertex position
+
+ TString GetCalorimeter() const {return fCalorimeter ; }
+ void SetCalorimeter(TString det) {fCalorimeter = det ; }
+
+ void Terminate();
+
+ private:
+ Bool_t IsBadRun(Int_t /*iRun*/) const {return kFALSE;} //Tests if this run bad according to private list
+
+ private:
+ Int_t fNCentrBin ; //number of bins in event container for centrality
+ Int_t fNZvertBin ; //number of bins in event container for vertex position
+ Int_t fNrpBin ; //number of bins in event container for reaction plain
+ Int_t fNPID ; //Number of possible PID combinations
+ Int_t fNmaxMixEv ; //Maximal number of events stored in buffer for mixing
+ Float_t fZvtxCut ; //Cut on vertex position
+ TString fCalorimeter ; //Select Calorimeter for IM
+
+ TList ** fEventsList ; //! containers for photons in stored events
+
+ //Histograms
+
+ TH3D * fhEtalon ; //Etalon histo, all distributions will have same binning as this one
+
+ TH3D ** fhRe1 ; //!REAL two-photon invariant mass distribution for different centralities and PID
+ TH3D ** fhMi1 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
+ TH3D ** fhRe2 ; //!REAL two-photon invariant mass distribution for different centralities and PID
+ TH3D ** fhMi2 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
+ TH3D ** fhRe3 ; //!REAL two-photon invariant mass distribution for different centralities and PID
+ TH3D ** fhMi3 ; //!MIXED two-photon invariant mass distribution for different centralities and PID
+ TH3D * fhEvents; //!Number of events per centrality, RP, zbin
+
+ //Acceptance
+ TH1D * fhPrimPt ; //! Spectrum of Primary
+ TH1D * fhPrimAccPt ; //! Spectrum of primary with accepted daughters
+ TH1D * fhPrimY ; //! Rapidity distribution of primary particles
+ TH1D * fhPrimAccY ; //! Rapidity distribution of primary with accepted daughters
+ TH1D * fhPrimPhi ; //! Azimutal distribution of primary particles
+ TH1D * fhPrimAccPhi; //! Azimutal distribution of primary with accepted daughters
+
#ifdef __PHOSGEO__
- AliPHOSGeoUtils * fPHOSGeo ; //! PHOS geometry pointer
+ AliPHOSGeoUtils * fPHOSGeo ; //! PHOS geometry pointer
#endif
-
- ClassDef(AliAnaPi0,3)
+
+ ClassDef(AliAnaPi0,3)
} ;
#include "AliIsolationCut.h"\r
#include "AliNeutralMesonSelection.h"\r
#include "AliCaloPID.h"\r
+#include "AliMCAnalysisUtils.h"\r
#include "AliAODPWG4ParticleCorrelation.h"\r
-\r
+ #include "AliStack.h"\r
ClassImp(AliAnaPi0EbE)\r
\r
//____________________________________________________________________________\r
AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));\r
mom1 = *(photon1->Momentum());\r
\r
- //Play with the MC stack if available\r
- \r
+ \r
for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){\r
\r
AliAODPWG4ParticleCorrelation * photon2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jphoton));\r
mom2 = *(photon2->Momentum());\r
+ \r
//Select good pair (good phi, pt cuts, aperture and invariant mass)\r
if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))\r
{\r
- if(GetDebug()>1) printf("Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
- \r
+ if(GetDebug()>1) \r
+ printf("Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
\r
+ //Play with the MC stack if available\r
if(IsDataMC()){\r
//Check origin of the candidates\r
- tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
- tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
+ tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
+ tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
+ \r
if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
- if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0;\r
+ if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){\r
+ \r
+ //Check if pi0 mother is the same\r
+ Int_t label1 = photon1->GetLabel();\r
+ TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree\r
+ label1 = mother1->GetFirstMother();\r
+ //mother1 = GetMCStack()->Particle(label1);//pi0\r
+ \r
+ Int_t label2 = photon2->GetLabel();\r
+ TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree\r
+ label2 = mother2->GetFirstMother();\r
+ //mother2 = GetMCStack()->Particle(label2);//pi0\r
+ \r
+ //printf("mother1 %d, mother2 %d\n",label1,label2);\r
+ if(label1 == label2)\r
+ tag = AliMCAnalysisUtils::kMCPi0;\r
+ }\r
}//Work with stack also \r
\r
//Create AOD for analysis\r
\r
if(IsDataMC()){\r
//Check origin of the candidates\r
- tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
- tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
+ tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
+ tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
- if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0;\r
+ if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){\r
+ //Check if pi0 mother is the same\r
+ Int_t label1 = photon1->GetLabel();\r
+ TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree\r
+ label1 = mother1->GetFirstMother();\r
+ //mother1 = GetMCStack()->Particle(label1);//pi0\r
+ \r
+ Int_t label2 = photon2->GetLabel();\r
+ TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree\r
+ label2 = mother2->GetFirstMother();\r
+ //mother2 = GetMCStack()->Particle(label2);//pi0\r
+ \r
+ //printf("mother1 %d, mother2 %d\n",label1,label2);\r
+ if(label1 == label2)\r
+ tag = AliMCAnalysisUtils::kMCPi0;\r
+ }\r
}//Work with stack also \r
\r
//Create AOD for analysis\r
//Play with the MC stack if available\r
//Check origin of the candidates\r
if(IsDataMC()){\r
- aodph.SetTag(GetCaloPID()->CheckOrigin(calo->GetLabel(0),GetMCStack()));\r
+ aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));\r
if(GetDebug() > 0) printf("AliAnaPi0EbE::FillAOD: Origin of candidate %d\n",aodph.GetTag());\r
}//Work with stack also \r
\r
fhEtaPi0 ->Fill(pt,eta);\r
\r
if(IsDataMC()){\r
- if(pi0->GetTag()== AliCaloPID::kMCPi0Decay){\r
+ if(pi0->GetTag()== AliMCAnalysisUtils::kMCPi0){\r
fhPtMCPi0 ->Fill(pt);\r
fhPhiMCPi0 ->Fill(pt,phi);\r
fhEtaMCPi0 ->Fill(pt,eta);\r
SRCS = PartCorrBase/AliAODPWG4Particle.cxx PartCorrBase/AliAODPWG4ParticleCorrelation.cxx \
PartCorrBase/AliNeutralMesonSelection.cxx PartCorrBase/AliFidutialCut.cxx \
- PartCorrBase/AliCaloPID.cxx PartCorrBase/AliIsolationCut.cxx\
+ PartCorrBase/AliCaloPID.cxx PartCorrBase/AliMCAnalysisUtils.cxx PartCorrBase/AliIsolationCut.cxx\
PartCorrBase/AliAnaScale.cxx PartCorrBase/AliAnaPartCorrMaker.cxx \
PartCorrBase/AliAnaPartCorrBaseClass.cxx PartCorrBase/AliAnalysisTaskParticleCorrelation.cxx \
PartCorrBase/AliCaloTrackReader.cxx PartCorrBase/AliCaloTrackESDReader.cxx \
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff