From 6639984f84fc5f87526d15c1597f1d54883765ff Mon Sep 17 00:00:00 2001 From: gconesab Date: Thu, 29 Jan 2009 10:28:56 +0000 Subject: [PATCH 1/1] 1)Terminate() method implemented in the frame. Simple examples on what to do with it put in AliAnaExample and AliAnaPi0 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 --- PWG4/CMake_libPWG4PartCorrBase.txt | 2 +- PWG4/PWG4PartCorrBaseLinkDef.h | 1 + PWG4/PartCorrBase/AliAnaPartCorrBaseClass.cxx | 24 +- PWG4/PartCorrBase/AliAnaPartCorrBaseClass.h | 17 +- PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx | 16 +- PWG4/PartCorrBase/AliAnaPartCorrMaker.h | 2 + .../AliAnalysisTaskParticleCorrelation.cxx | 6 +- PWG4/PartCorrBase/AliCaloPID.cxx | 118 +----- PWG4/PartCorrBase/AliCaloPID.h | 9 +- PWG4/PartCorrBase/AliCaloTrackMCReader.cxx | 13 +- PWG4/PartCorrBase/AliCaloTrackMCReader.h | 2 +- PWG4/PartCorrBase/AliCaloTrackReader.cxx | 47 +-- PWG4/PartCorrBase/AliCaloTrackReader.h | 7 +- PWG4/PartCorrBase/AliMCAnalysisUtils.cxx | 355 ++++++++++++++++++ PWG4/PartCorrBase/AliMCAnalysisUtils.h | 63 ++++ PWG4/PartCorrDep/AliAnaCaloTrigger.cxx | 281 +++++++------- PWG4/PartCorrDep/AliAnaCaloTrigger.h | 24 +- PWG4/PartCorrDep/AliAnaExample.cxx | 44 ++- PWG4/PartCorrDep/AliAnaExample.h | 2 + .../AliAnaParticleHadronCorrelation.cxx | 15 +- .../AliAnaParticleHadronCorrelation.h | 8 +- PWG4/PartCorrDep/AliAnaParticleIsolation.cxx | 42 +-- PWG4/PartCorrDep/AliAnaPhoton.cxx | 97 +++-- PWG4/PartCorrDep/AliAnaPhoton.h | 9 +- PWG4/PartCorrDep/AliAnaPi0.cxx | 95 ++++- PWG4/PartCorrDep/AliAnaPi0.h | 149 ++++---- PWG4/PartCorrDep/AliAnaPi0EbE.cxx | 60 ++- PWG4/libPWG4PartCorrBase.pkg | 2 +- 28 files changed, 1040 insertions(+), 470 deletions(-) create mode 100755 PWG4/PartCorrBase/AliMCAnalysisUtils.cxx create mode 100755 PWG4/PartCorrBase/AliMCAnalysisUtils.h diff --git a/PWG4/CMake_libPWG4PartCorrBase.txt b/PWG4/CMake_libPWG4PartCorrBase.txt index 4fd80f3181d..91e09924dbe 100755 --- a/PWG4/CMake_libPWG4PartCorrBase.txt +++ b/PWG4/CMake_libPWG4PartCorrBase.txt @@ -3,7 +3,7 @@ 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 diff --git a/PWG4/PWG4PartCorrBaseLinkDef.h b/PWG4/PWG4PartCorrBaseLinkDef.h index 316ae4cd2e9..cd71e0ed5cc 100755 --- a/PWG4/PWG4PartCorrBaseLinkDef.h +++ b/PWG4/PWG4PartCorrBaseLinkDef.h @@ -10,6 +10,7 @@ #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+; diff --git a/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.cxx b/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.cxx index ed52d4aae64..c4ea5c239ea 100755 --- a/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.cxx +++ b/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.cxx @@ -33,6 +33,7 @@ #include "AliCaloPID.h" #include "AliFidutialCut.h" #include "AliIsolationCut.h" +#include "AliMCAnalysisUtils.h" #include "AliNeutralMesonSelection.h" #include "AliLog.h" #include "AliAODPWG4ParticleCorrelation.h" @@ -60,6 +61,7 @@ ClassImp(AliAnaPartCorrBaseClass) fCaloPID = new AliCaloPID(); fFidCut = new AliFidutialCut(); fIC = new AliIsolationCut(); + fMCUtils = new AliMCAnalysisUtils(); //Initialize parameters InitParameters(); @@ -76,7 +78,7 @@ AliAnaPartCorrBaseClass::AliAnaPartCorrBaseClass(const AliAnaPartCorrBaseClass & 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) @@ -108,8 +110,9 @@ AliAnaPartCorrBaseClass & AliAnaPartCorrBaseClass::operator = (const AliAnaPartC 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) ; @@ -117,7 +120,6 @@ AliAnaPartCorrBaseClass & AliAnaPartCorrBaseClass::operator = (const AliAnaPartC 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; @@ -151,11 +153,12 @@ AliAnaPartCorrBaseClass::~AliAnaPartCorrBaseClass() 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 ; } @@ -303,6 +306,13 @@ TNamed * AliAnaPartCorrBaseClass::GetEMCALCells() const { } +//__________________________________________________ +Int_t AliAnaPartCorrBaseClass::GetEventNumber() const { + //Get current event number + + return fReader->GetEventNumber() ; +} + //__________________________________________________ AliStack * AliAnaPartCorrBaseClass::GetMCStack() const { //Get stack pointer from reader diff --git a/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.h b/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.h index 3766480ad30..3dbdb350f2b 100755 --- a/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.h +++ b/PWG4/PartCorrBase/AliAnaPartCorrBaseClass.h @@ -21,6 +21,7 @@ class AliCaloTrackReader ; class AliCaloPID ; class AliFidutialCut ; class AliIsolationCut ; +class AliMCAnalysisUtils ; class AliNeutralMesonSelection ; class AliStack ; class AliHeader ; @@ -58,12 +59,15 @@ public: 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 ; } @@ -91,7 +95,7 @@ public: 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 ;} @@ -101,6 +105,9 @@ public: 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 ;} @@ -127,7 +134,6 @@ public: 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 ; @@ -181,9 +187,12 @@ public: 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 @@ -197,7 +206,7 @@ public: Float_t fHistoEtaMax ; //Maximum value of eta histogram range Float_t fHistoEtaMin ; //Minimum value of eta histogram range - ClassDef(AliAnaPartCorrBaseClass,2) + ClassDef(AliAnaPartCorrBaseClass,3) } ; diff --git a/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx b/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx index 22104e898dd..17f1cf50d7b 100755 --- a/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx +++ b/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx @@ -198,7 +198,7 @@ Bool_t AliAnaPartCorrMaker::ProcessEvent(Int_t iEntry){ 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"); @@ -223,3 +223,17 @@ Bool_t AliAnaPartCorrMaker::ProcessEvent(Int_t iEntry){ 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 +} diff --git a/PWG4/PartCorrBase/AliAnaPartCorrMaker.h b/PWG4/PartCorrBase/AliAnaPartCorrMaker.h index f453e78f17c..c7bf25c8091 100755 --- a/PWG4/PartCorrBase/AliAnaPartCorrMaker.h +++ b/PWG4/PartCorrBase/AliAnaPartCorrMaker.h @@ -46,6 +46,8 @@ public: 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"); diff --git a/PWG4/PartCorrBase/AliAnalysisTaskParticleCorrelation.cxx b/PWG4/PartCorrBase/AliAnalysisTaskParticleCorrelation.cxx index c3a0ceacf0d..0a7172f885c 100755 --- a/PWG4/PartCorrBase/AliAnalysisTaskParticleCorrelation.cxx +++ b/PWG4/PartCorrBase/AliAnalysisTaskParticleCorrelation.cxx @@ -171,8 +171,10 @@ void AliAnalysisTaskParticleCorrelation::UserExec(Option_t */*option*/) void AliAnalysisTaskParticleCorrelation::Terminate(Option_t */*option*/) { // Terminate analysis + // Do some plots // - AliDebug(1,"Do nothing in Terminate"); - //fAna->Terminate(); + + fAna->Terminate(); + } diff --git a/PWG4/PartCorrBase/AliCaloPID.cxx b/PWG4/PartCorrBase/AliCaloPID.cxx index 023f0399b45..ee0efe8c47c 100755 --- a/PWG4/PartCorrBase/AliCaloPID.cxx +++ b/PWG4/PartCorrBase/AliCaloPID.cxx @@ -33,8 +33,6 @@ #include "AliCaloPID.h" #include "AliAODCaloCluster.h" #include "AliAODPWG4Particle.h" -#include "AliStack.h" -#include "TParticle.h" ClassImp(AliCaloPID) @@ -48,7 +46,7 @@ fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.), 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 @@ -72,7 +70,7 @@ fPHOSWeightFormula(pid.fPHOSWeightFormula), fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula), fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula), fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut), -fDebug(pid.fDebug),fMCGenerator(pid.fMCGenerator) +fDebug(pid.fDebug) { // cpy ctor @@ -104,8 +102,7 @@ AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid) fDispCut = pid.fDispCut; fTOFCut = pid.fTOFCut; fDebug = pid.fDebug; - fMCGenerator = pid.fMCGenerator; - + return *this; } @@ -120,99 +117,6 @@ AliCaloPID::~AliCaloPID() { } -//_________________________________________________________________________ -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() { @@ -240,7 +144,6 @@ void AliCaloPID::InitParameters() fDispCut = 1.5; fTOFCut = 5.e-9; fDebug = -1; - fMCGenerator = "PYTHIA"; } //_______________________________________________________________ @@ -405,7 +308,6 @@ void AliCaloPID::Print(const Option_t * opt) const 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"); @@ -431,13 +333,13 @@ void AliCaloPID::SetPIDBits(const TString calo, const AliAODCaloCluster * cluste 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()); + } } diff --git a/PWG4/PartCorrBase/AliCaloPID.h b/PWG4/PartCorrBase/AliCaloPID.h index 519069918e4..9c36bb65b3b 100755 --- a/PWG4/PartCorrBase/AliCaloPID.h +++ b/PWG4/PartCorrBase/AliCaloPID.h @@ -20,7 +20,6 @@ class TFormula ; class AliLog ; class AliAODCaloCluster; class AliAODPWG4Particle; -#include "AliStack.h" class AliCaloPID : public TObject { @@ -43,11 +42,9 @@ public: 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 ; @@ -101,9 +98,6 @@ public: 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 @@ -125,7 +119,6 @@ private: 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) } ; diff --git a/PWG4/PartCorrBase/AliCaloTrackMCReader.cxx b/PWG4/PartCorrBase/AliCaloTrackMCReader.cxx index 85d44282103..05c9c2df13f 100755 --- a/PWG4/PartCorrBase/AliCaloTrackMCReader.cxx +++ b/PWG4/PartCorrBase/AliCaloTrackMCReader.cxx @@ -173,9 +173,10 @@ void AliCaloTrackMCReader::FillCalorimeters(const Int_t iParticle, TParticle* p } //____________________________________________________________________________ -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); @@ -190,14 +191,14 @@ void AliCaloTrackMCReader::FillInputEvent() 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]; diff --git a/PWG4/PartCorrBase/AliCaloTrackMCReader.h b/PWG4/PartCorrBase/AliCaloTrackMCReader.h index 48077d69631..ab879b8c5c1 100755 --- a/PWG4/PartCorrBase/AliCaloTrackMCReader.h +++ b/PWG4/PartCorrBase/AliCaloTrackMCReader.h @@ -59,7 +59,7 @@ class AliCaloTrackMCReader : public AliCaloTrackReader { 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) ; diff --git a/PWG4/PartCorrBase/AliCaloTrackReader.cxx b/PWG4/PartCorrBase/AliCaloTrackReader.cxx index 1e85be8e641..a61d8caa5a2 100755 --- a/PWG4/PartCorrBase/AliCaloTrackReader.cxx +++ b/PWG4/PartCorrBase/AliCaloTrackReader.cxx @@ -43,7 +43,8 @@ ClassImp(AliCaloTrackReader) //____________________________________________________________________________ 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), @@ -59,7 +60,8 @@ ClassImp(AliCaloTrackReader) //____________________________________________________________________________ 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)), @@ -81,30 +83,30 @@ AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & 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; @@ -163,7 +165,7 @@ AliHeader* AliCaloTrackReader::GetHeader() const { if(fMC) return fMC->Header(); else{ - printf("header is not available"); + printf("Header is not available\n"); return 0x0 ; } } @@ -173,7 +175,7 @@ AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const { if(fMC) return fMC->GenEventHeader(); else{ - printf("GenEventHeader is not available"); + printf("GenEventHeader is not available\n"); return 0x0 ; } } @@ -222,9 +224,10 @@ void AliCaloTrackReader::Print(const Option_t * opt) const } //___________________________________________________ -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(); diff --git a/PWG4/PartCorrBase/AliCaloTrackReader.h b/PWG4/PartCorrBase/AliCaloTrackReader.h index 7d383759a75..5bcf7367980 100755 --- a/PWG4/PartCorrBase/AliCaloTrackReader.h +++ b/PWG4/PartCorrBase/AliCaloTrackReader.h @@ -50,6 +50,7 @@ public: 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 ; } @@ -82,7 +83,7 @@ public: 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() {;} @@ -116,11 +117,11 @@ public: 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 diff --git a/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx b/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx new file mode 100755 index 00000000000..e0dcff2e5c3 --- /dev/null +++ b/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx @@ -0,0 +1,355 @@ +/************************************************************************** + * 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 +#include + +//---- 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"); + +} + + diff --git a/PWG4/PartCorrBase/AliMCAnalysisUtils.h b/PWG4/PartCorrBase/AliMCAnalysisUtils.h new file mode 100755 index 00000000000..a6484227f70 --- /dev/null +++ b/PWG4/PartCorrBase/AliMCAnalysisUtils.h @@ -0,0 +1,63 @@ +#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 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 + + + diff --git a/PWG4/PartCorrDep/AliAnaCaloTrigger.cxx b/PWG4/PartCorrDep/AliAnaCaloTrigger.cxx index 140ce55cd51..a941a8d9ed2 100755 --- a/PWG4/PartCorrDep/AliAnaCaloTrigger.cxx +++ b/PWG4/PartCorrDep/AliAnaCaloTrigger.cxx @@ -18,25 +18,24 @@ // 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 -#include +//Root #include #include +//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), @@ -49,9 +48,7 @@ AliAnaCaloTrigger::AliAnaCaloTrigger() : //______________________________________________________________________________ AliAnaCaloTrigger::AliAnaCaloTrigger(const char *name) : - AliAnalysisTask(name,"AnaCaloTrigger"), - fChain(0), - fESD(0), + AliAnalysisTaskSE(name), fOutputContainer(0), fCalorimeter("PHOS"), fNtTrigger22(0), @@ -59,22 +56,19 @@ AliAnaCaloTrigger::AliAnaCaloTrigger(const char *name) : { // 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()) ; - + } //_________________________________________________________________________ @@ -82,10 +76,10 @@ AliAnaCaloTrigger & AliAnaCaloTrigger::operator = (const AliAnaCaloTrigger & sou { // 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 ; @@ -99,46 +93,28 @@ AliAnaCaloTrigger & AliAnaCaloTrigger::operator = (const AliAnaCaloTrigger & sou 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(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) ; @@ -147,107 +123,118 @@ void AliAnaCaloTrigger::CreateOutputObjects() } //______________________________________________________________________________ -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(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(triggerAmplitudes->At(0)) ; - const Float_t ka22O = static_cast(triggerAmplitudes->At(1)) ; - const Float_t kaNN = static_cast(triggerAmplitudes->At(2)) ; - const Float_t kaNNO = static_cast(triggerAmplitudes->At(3)) ; - - // trigger position - const Float_t kx22 = static_cast(triggerPosition->At(0)) ; - const Float_t ky22 = static_cast(triggerPosition->At(1)) ; - const Float_t kz22 = static_cast(triggerPosition->At(2)) ; - const Float_t kxNN = static_cast(triggerPosition->At(3)) ; - const Float_t kyNN = static_cast(triggerPosition->At(4)) ; - const Float_t kzNN = static_cast(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(triggerAmplitudes->At(0)) ; + const Float_t ka22O = static_cast(triggerAmplitudes->At(1)) ; + const Float_t kaNN = static_cast(triggerAmplitudes->At(2)) ; + const Float_t kaNNO = static_cast(triggerAmplitudes->At(3)) ; + + // trigger position + const Float_t kx22 = static_cast(triggerPosition->At(0)) ; + const Float_t ky22 = static_cast(triggerPosition->At(1)) ; + const Float_t kz22 = static_cast(triggerPosition->At(2)) ; + const Float_t kxNN = static_cast(triggerPosition->At(3)) ; + const Float_t kyNN = static_cast(triggerPosition->At(4)) ; + const Float_t kzNN = static_cast(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); + } //______________________________________________________________________________ diff --git a/PWG4/PartCorrDep/AliAnaCaloTrigger.h b/PWG4/PartCorrDep/AliAnaCaloTrigger.h index c173eacbaff..7c1a660c7fa 100755 --- a/PWG4/PartCorrDep/AliAnaCaloTrigger.h +++ b/PWG4/PartCorrDep/AliAnaCaloTrigger.h @@ -7,22 +7,18 @@ // 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() ; @@ -31,26 +27,22 @@ public: 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 diff --git a/PWG4/PartCorrDep/AliAnaExample.cxx b/PWG4/PartCorrDep/AliAnaExample.cxx index f77a62eff82..b26c5ea1f67 100755 --- a/PWG4/PartCorrDep/AliAnaExample.cxx +++ b/PWG4/PartCorrDep/AliAnaExample.cxx @@ -30,6 +30,9 @@ //#include "Riostream.h" #include "TClonesArray.h" #include "TParticle.h" +#include "TCanvas.h" +#include "TPad.h" +#include "TROOT.h" //---- AliRoot system ---- #include "AliAnaExample.h" @@ -316,8 +319,8 @@ void AliAnaExample::MakeAnalysisFillAOD() 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()); } } @@ -372,3 +375,40 @@ void AliAnaExample::MakeAnalysisFillHistograms() }//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()); + +} diff --git a/PWG4/PartCorrDep/AliAnaExample.h b/PWG4/PartCorrDep/AliAnaExample.h index cac88c15d37..a034079c6b0 100755 --- a/PWG4/PartCorrDep/AliAnaExample.h +++ b/PWG4/PartCorrDep/AliAnaExample.h @@ -42,6 +42,8 @@ class AliAnaExample : public AliAnaPartCorrBaseClass { TString GetDetector() const {return fDetector ;} void SetDetector( TString calo ) {fDetector = calo; } + void Terminate(); + private: Int_t fPdg ; //identified particle id diff --git a/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.cxx b/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.cxx index bc5c925ea37..1e8607ab6b0 100755 --- a/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.cxx +++ b/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.cxx @@ -38,8 +38,8 @@ ClassImp(AliAnaParticleHadronCorrelation) //____________________________________________________________________________ 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), @@ -55,6 +55,7 @@ ClassImp(AliAnaParticleHadronCorrelation) 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), @@ -80,7 +81,8 @@ AliAnaParticleHadronCorrelation & AliAnaParticleHadronCorrelation::operator = (c fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ; fDeltaPhiMinCut = source.fDeltaPhiMinCut ; - + fSelectIsolated = source.fSelectIsolated ; + fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ; fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ; fhDeltaPhiCharged = source.fhDeltaPhiCharged ; @@ -230,7 +232,7 @@ void AliAnaParticleHadronCorrelation::InitParameters() SetPtCutRange(2,300); fDeltaPhiMinCut = 1.5 ; fDeltaPhiMaxCut = 4.5 ; - + fSelectIsolated = kFALSE; } //__________________________________________________________________ @@ -246,8 +248,7 @@ void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const 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) ; } //____________________________________________________________________________ @@ -318,6 +319,8 @@ void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() 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); diff --git a/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.h b/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.h index 27333aa4604..dab3d945a92 100755 --- a/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.h +++ b/PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.h @@ -32,6 +32,9 @@ public: 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; @@ -48,8 +51,9 @@ public: 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 diff --git a/PWG4/PartCorrDep/AliAnaParticleIsolation.cxx b/PWG4/PartCorrDep/AliAnaParticleIsolation.cxx index 293539d4a07..3e937102b32 100755 --- a/PWG4/PartCorrDep/AliAnaParticleIsolation.cxx +++ b/PWG4/PartCorrDep/AliAnaParticleIsolation.cxx @@ -35,7 +35,7 @@ #include "AliIsolationCut.h" #include "AliNeutralMesonSelection.h" #include "AliAODPWG4ParticleCorrelation.h" -#include "AliCaloPID.h" +#include "AliMCAnalysisUtils.h" ClassImp(AliAnaParticleIsolation) @@ -792,30 +792,30 @@ void AliAnaParticleIsolation::MakeAnalysisFillHistograms() if(IsDataMC()){ Int_t tag =aod->GetTag(); - if(tag == AliCaloPID::kMCPrompt){ + if(tag == AliMCAnalysisUtils::kMCPrompt){ fhPtIsoPrompt ->Fill(ptcluster); fhPhiIsoPrompt ->Fill(ptcluster,phicluster); fhEtaIsoPrompt ->Fill(ptcluster,etacluster); } - else if(tag == AliCaloPID::kMCFragmentation) + else if(tag == AliMCAnalysisUtils::kMCFragmentation) { fhPtIsoFragmentation ->Fill(ptcluster); fhPhiIsoFragmentation ->Fill(ptcluster,phicluster); fhEtaIsoFragmentation ->Fill(ptcluster,etacluster); } - else if(tag == AliCaloPID::kMCPi0Decay) + else if(tag == AliMCAnalysisUtils::kMCPi0Decay) { fhPtIsoPi0Decay ->Fill(ptcluster); fhPhiIsoPi0Decay ->Fill(ptcluster,phicluster); fhEtaIsoPi0Decay ->Fill(ptcluster,etacluster); } - else if(tag == AliCaloPID::kMCEtaDecay || tag == AliCaloPID::kMCOtherDecay) + else if(tag == AliMCAnalysisUtils::kMCEtaDecay || tag == AliMCAnalysisUtils::kMCOtherDecay) { fhPtIsoOtherDecay ->Fill(ptcluster); fhPhiIsoOtherDecay ->Fill(ptcluster,phicluster); fhEtaIsoOtherDecay ->Fill(ptcluster,etacluster); } - else if(tag == AliCaloPID::kMCConversion) + else if(tag == AliMCAnalysisUtils::kMCConversion) { fhPtIsoConversion ->Fill(ptcluster); fhPhiIsoConversion ->Fill(ptcluster,phicluster); @@ -899,11 +899,11 @@ void AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati if(n[icone][ipt] == 0) { fhPtThresIsolated[icone][ipt]->Fill(ptC); if(IsDataMC()){ - if(tag == AliCaloPID::kMCPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCOtherDecay || tag == AliCaloPID::kMCEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC) ; + if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC) ; else fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC) ; } } @@ -912,11 +912,11 @@ void AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati if(nfrac[icone][ipt] == 0) { fhPtFracIsolated[icone][ipt]->Fill(ptC); if(IsDataMC()){ - if(tag == AliCaloPID::kMCPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ; - else if(tag == AliCaloPID::kMCOtherDecay || tag == AliCaloPID::kMCEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC) ; + if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC) ; + else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC) ; else fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC) ; } } @@ -925,11 +925,11 @@ void AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati //Sum in cone histograms fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ; if(IsDataMC()){ - if(tag == AliCaloPID::kMCPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ; - else if(tag == AliCaloPID::kMCConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ; - else if(tag == AliCaloPID::kMCFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ; - else if(tag == AliCaloPID::kMCPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ; - else if(tag == AliCaloPID::kMCOtherDecay || tag == AliCaloPID::kMCEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ; + if(tag == AliMCAnalysisUtils::kMCPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ; + else if(tag == AliMCAnalysisUtils::kMCConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ; + else if(tag == AliMCAnalysisUtils::kMCFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ; + else if(tag == AliMCAnalysisUtils::kMCPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ; + else if(tag == AliMCAnalysisUtils::kMCOtherDecay || tag == AliMCAnalysisUtils::kMCEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ; else fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ; } diff --git a/PWG4/PartCorrDep/AliAnaPhoton.cxx b/PWG4/PartCorrDep/AliAnaPhoton.cxx index 10327487b1a..655c6e1faf2 100755 --- a/PWG4/PartCorrDep/AliAnaPhoton.cxx +++ b/PWG4/PartCorrDep/AliAnaPhoton.cxx @@ -31,6 +31,8 @@ #include "AliAnaPhoton.h" #include "AliCaloTrackReader.h" #include "AliCaloPID.h" +#include "AliMCAnalysisUtils.h" +#include "AliStack.h" ClassImp(AliAnaPhoton) @@ -42,6 +44,7 @@ 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), @@ -62,10 +65,12 @@ AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) : //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 @@ -83,26 +88,29 @@ AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g) 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; @@ -172,7 +180,7 @@ TList * AliAnaPhoton::GetCreateOutputObjects() 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)"); @@ -190,6 +198,23 @@ TList * AliAnaPhoton::GetCreateOutputObjects() 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)"); @@ -240,7 +265,7 @@ TList * AliAnaPhoton::GetCreateOutputObjects() 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)"); @@ -257,6 +282,7 @@ TList * AliAnaPhoton::GetCreateOutputObjects() fhEtaUnknown->SetYTitle("#eta"); fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)"); outputContainer->Add(fhEtaUnknown) ; + }//Histos with MC //Save parameters used for analysis @@ -290,6 +316,24 @@ TList * AliAnaPhoton::GetCreateOutputObjects() } +//____________________________________________________________________________ +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() { @@ -397,7 +441,7 @@ void AliAnaPhoton::MakeAnalysisFillAOD() //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 @@ -447,40 +491,47 @@ void AliAnaPhoton::MakeAnalysisFillHistograms() 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); } diff --git a/PWG4/PartCorrDep/AliAnaPhoton.h b/PWG4/PartCorrDep/AliAnaPhoton.h index d913dba2274..ebe3af71872 100755 --- a/PWG4/PartCorrDep/AliAnaPhoton.h +++ b/PWG4/PartCorrDep/AliAnaPhoton.h @@ -33,6 +33,8 @@ public: TList * GetCreateOutputObjects(); + void Init(); + void MakeAnalysisFillAOD() ; void MakeAnalysisFillHistograms() ; @@ -71,6 +73,10 @@ public: 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 @@ -87,7 +93,8 @@ public: TH2F * fhPhiUnknown; //! Phi of identified Unknown gamma TH2F * fhEtaUnknown; //! eta of identified Unknown gamma - ClassDef(AliAnaPhoton,1) + ClassDef(AliAnaPhoton,2) + } ; diff --git a/PWG4/PartCorrDep/AliAnaPi0.cxx b/PWG4/PartCorrDep/AliAnaPi0.cxx index a7195a68f13..18537651b73 100755 --- a/PWG4/PartCorrDep/AliAnaPi0.cxx +++ b/PWG4/PartCorrDep/AliAnaPi0.cxx @@ -27,11 +27,15 @@ // --- 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" @@ -304,7 +308,7 @@ void AliAnaPi0::Print(const Option_t * /*opt*/) const //____________________________________________________________________________________________________________________________________________________ -void AliAnaPi0:: MakeAnalysisFillHistograms() +void AliAnaPi0::MakeAnalysisFillHistograms() { //Process one event and extract photons from AOD branch // filled with AliAnaPhoton and fill histos with invariant mass @@ -471,8 +475,97 @@ void AliAnaPi0:: MakeAnalysisFillHistograms() +//____________________________________________________________________________________________________________________________________________________ +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()); +} diff --git a/PWG4/PartCorrDep/AliAnaPi0.h b/PWG4/PartCorrDep/AliAnaPi0.h index 3c7cdfb7992..dbb18c28ed0 100755 --- a/PWG4/PartCorrDep/AliAnaPi0.h +++ b/PWG4/PartCorrDep/AliAnaPi0.h @@ -22,85 +22,86 @@ class AliESDEvent ; #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) } ; diff --git a/PWG4/PartCorrDep/AliAnaPi0EbE.cxx b/PWG4/PartCorrDep/AliAnaPi0EbE.cxx index b8ed86f4492..67d4b59780f 100755 --- a/PWG4/PartCorrDep/AliAnaPi0EbE.cxx +++ b/PWG4/PartCorrDep/AliAnaPi0EbE.cxx @@ -36,8 +36,9 @@ #include "AliIsolationCut.h" #include "AliNeutralMesonSelection.h" #include "AliCaloPID.h" +#include "AliMCAnalysisUtils.h" #include "AliAODPWG4ParticleCorrelation.h" - + #include "AliStack.h" ClassImp(AliAnaPi0EbE) //____________________________________________________________________________ @@ -264,24 +265,42 @@ void AliAnaPi0EbE::MakeInvMassInCalorimeter() AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton)); mom1 = *(photon1->Momentum()); - //Play with the MC stack if available - + for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){ AliAODPWG4ParticleCorrelation * photon2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jphoton)); mom2 = *(photon2->Momentum()); + //Select good pair (good phi, pt cuts, aperture and invariant mass) if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)) { - 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()); - + 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()); + //Play with the MC stack if available if(IsDataMC()){ //Check origin of the candidates - tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack()); - tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack()); + tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack()); + tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack()); + if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2); - if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0; + if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){ + + //Check if pi0 mother is the same + Int_t label1 = photon1->GetLabel(); + TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree + label1 = mother1->GetFirstMother(); + //mother1 = GetMCStack()->Particle(label1);//pi0 + + Int_t label2 = photon2->GetLabel(); + TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree + label2 = mother2->GetFirstMother(); + //mother2 = GetMCStack()->Particle(label2);//pi0 + + //printf("mother1 %d, mother2 %d\n",label1,label2); + if(label1 == label2) + tag = AliMCAnalysisUtils::kMCPi0; + } }//Work with stack also //Create AOD for analysis @@ -339,10 +358,25 @@ void AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() if(IsDataMC()){ //Check origin of the candidates - tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack()); - tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack()); + tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack()); + tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack()); if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2); - if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0; + if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){ + //Check if pi0 mother is the same + Int_t label1 = photon1->GetLabel(); + TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree + label1 = mother1->GetFirstMother(); + //mother1 = GetMCStack()->Particle(label1);//pi0 + + Int_t label2 = photon2->GetLabel(); + TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree + label2 = mother2->GetFirstMother(); + //mother2 = GetMCStack()->Particle(label2);//pi0 + + //printf("mother1 %d, mother2 %d\n",label1,label2); + if(label1 == label2) + tag = AliMCAnalysisUtils::kMCPi0; + } }//Work with stack also //Create AOD for analysis @@ -454,7 +488,7 @@ void AliAnaPi0EbE::MakeShowerShapeIdentification() //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("AliAnaPi0EbE::FillAOD: Origin of candidate %d\n",aodph.GetTag()); }//Work with stack also @@ -496,7 +530,7 @@ void AliAnaPi0EbE::MakeAnalysisFillHistograms() fhEtaPi0 ->Fill(pt,eta); if(IsDataMC()){ - if(pi0->GetTag()== AliCaloPID::kMCPi0Decay){ + if(pi0->GetTag()== AliMCAnalysisUtils::kMCPi0){ fhPtMCPi0 ->Fill(pt); fhPhiMCPi0 ->Fill(pt,phi); fhEtaMCPi0 ->Fill(pt,eta); diff --git a/PWG4/libPWG4PartCorrBase.pkg b/PWG4/libPWG4PartCorrBase.pkg index 030e7f5c72c..06f5e3827cb 100755 --- a/PWG4/libPWG4PartCorrBase.pkg +++ b/PWG4/libPWG4PartCorrBase.pkg @@ -2,7 +2,7 @@ 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 \ -- 2.43.0