#include "TList.h"
#include "TH1F.h"
#include "TH2F.h"
+#include "TH1I.h"
#include "TTree.h"
#include <iostream>
#include "AliAnalysisEtCuts.h"
#include "Rtypes.h"
#include "TString.h"
#include "AliCentrality.h"
+#include "AliAnalysisEtSelector.h"
+
//#include "THnSparse.h"
using namespace std;
,fTrackDistanceCut(0)
,fTrackDxCut(0)
,fTrackDzCut(0)
+ ,fChargedEnergyRemoved(0)
+ ,fNeutralEnergyRemoved(0)
+ ,fGammaEnergyAdded(0)
,fHistEt(0)
,fHistChargedEt(0)
,fHistNeutralEt(0)
,fTreeDeposit(0)
,fCentrality(0)
,fDetector(0)
- ,fMakeSparse(kFALSE)
+ ,fMakeSparse(kTRUE)
,fSparseHistTracks(0)
,fSparseHistClusters(0)
,fSparseHistEt(0)
- ,fSparseTracks(0)
+ ,fSparseTracks(0)
,fSparseClusters(0)
,fSparseEt(0)
+ ,fClusterType(0)
+ ,fMatrixInitialized(kFALSE)
+ ,fCutFlow(0)
+ ,fSelector(0)
{}
list->Add(fSparseHistClusters);
list->Add(fSparseHistEt);
}
+
+ list->Add(fCutFlow);
}
fSparseHistEt->GetAxis(5)->SetTitle("charged_mult");
fSparseHistEt->GetAxis(6)->SetTitle("cent");
}
+
+ histname = "fCutFlow" + fHistogramNameSuffix;
+ fCutFlow = new TH1I(histname, histname, 20, -0.5, 19.5);
}
TH2F* AliAnalysisEt::CreateEtaEHisto2D(TString name, TString title, TString ztitle)
fTree->Branch("fChargedMultiplicity",&fChargedMultiplicity,"fChargedMultiplicity/I");
fTree->Branch("fNeutralMultiplicity",&fNeutralMultiplicity,"fNeutralMultiplicity/I");
fTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
+ fTree->Branch("fChargedEnergyRemoved", &fChargedEnergyRemoved, "fChargedEnergyRemoved/D");
+ fTree->Branch("fNeutralEnergyRemoved", &fNeutralEnergyRemoved, "fNeutralEnergyRemoved/D");
+ fTree->Branch("fGammaEnergyAdded", &fGammaEnergyAdded, "fGammaEnergyAdded/D");
fTree->Branch("fBaryonEt",&fBaryonEt,"fBaryonEt/D");
fTree->Branch("fAntiBaryonEt",&fAntiBaryonEt,"fAntiBaryonEt/D");
Int_t AliAnalysisEt::AnalyseEvent(AliVEvent *event)
{ //this line is basically here to eliminate a compiler warning that event is not used. Making it a virtual function did not work with the plugin.
+ fSelector->SetEvent(event);
AliAnalysisEtCommon::AnalyseEvent(event);
ResetEventValues();
return 0;
Float_t pos[3];
cluster->GetPosition(pos);
TVector3 cp(pos);
+ Double_t corrEnergy = 0;
+
+ if(cluster->E() < 1.5)
+ {
+ corrEnergy =cluster->E()/(0.51 + 0.02*cluster->E());
+ }
+ else
+ {
+ corrEnergy =cluster->E()/(0.51 + 0.02*1.5);
+ }
+ //std::cout << "Original energy: " << cluster->E() << ", corrected energy: " << corrEnergy << std::endl;
- return cluster->E() * TMath::Sin(cp.Theta());
+ return corrEnergy * TMath::Sin(cp.Theta());
+}
+
+//____________________________________________________________________________
+Double_t AliAnalysisEt::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, AliVEvent* e) const {
+ //Parameterization of LHC10h period
+ //_true if neutral_
+
+ Double_t meanX=0;
+ Double_t meanZ=0.;
+ Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
+ 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
+ Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
+
+ Double_t mf = e->GetMagneticField(); //Positive for ++ and negative for --
+
+ if(mf<0.){ //field --
+ meanZ = -0.468318 ;
+ if(charge>0)
+ meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
+ else
+ meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
+ }
+ else{ //Field ++
+ meanZ= -0.468318;
+ if(charge>0)
+ meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
+ else
+ meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
+ }
+ Double_t rz=(dz-meanZ)/sz ;
+ Double_t rx=(dx-meanX)/sx ;
+ return TMath::Sqrt(rx*rx+rz*rz) ;
}
#ifndef ALIANALYSISET_H
#define ALIANALYSISET_H
-class AliCentrality;
#include "AliAnalysisEtCommon.h"
#include "THnSparse.h"
+#include "AliESDCaloCluster.h"
+#include "AliAnalysisEtCuts.h"
+#include <vector>
+#include "Rtypes.h"
+class AliAnalysisEtSelector;
+class AliCentrality;
class TString;
class TTree;
class TH2F;
class TH1F;
+class TH1I;
class AliVEvent;
class TList;
class TString;
class AliAnalysisEt : public AliAnalysisEtCommon
{
public:
-
+
AliAnalysisEt();
virtual ~AliAnalysisEt();
-
+
public:
-
+
/** Analyse the event! */
virtual Int_t AnalyseEvent(AliVEvent *event);
/** Initialise the analysis, must be overloaded. */
virtual void Init();
- /**
- * Creates the histograms, must be overloaded if you want to add your own.
+ /**
+ * Creates the histograms, must be overloaded if you want to add your own.
* Uses the fHistogramNameSuffix to create proper histogram names
*/
virtual void CreateHistograms();
virtual void CreateTrees();
TH2F* CreateEtaEHisto2D(TString name, TString title, TString ztitle);
- TH2F* CreateEtaPtHisto2D(TString name, TString title, TString ztitle);
- TH2F* CreateEtaEtHisto2D(TString name, TString title, TString ztitle);
- TH2F* CreateResEHisto2D(TString name, TString title, TString ztitle);
- TH2F* CreateResPtHisto2D(TString name, TString title, TString ztitle);
+ TH2F* CreateEtaPtHisto2D(TString name, TString title, TString ztitle);
+ TH2F* CreateEtaEtHisto2D(TString name, TString title, TString ztitle);
+ TH2F* CreateResEHisto2D(TString name, TString title, TString ztitle);
+ TH2F* CreateResPtHisto2D(TString name, TString title, TString ztitle);
THnSparseF* CreateClusterHistoSparse(TString name, TString title);
THnSparseF* CreateNeutralPartHistoSparse(TString name, TString title);
THnSparseF* CreateChargedPartHistoSparse(TString name, TString title);
-
+
/** Fills the histograms, must be overloaded if you want to add your own */
virtual void FillHistograms();
virtual void ResetEventValues();
/** Total Et in the event (without acceptance cuts) */
- Double_t GetTotEt() const { return fTotEt; }
+ Double_t GetTotEt() const {
+ return fTotEt;
+ }
/** Total Et in the event within the acceptance cuts */
- Double_t GetTotEtAcc() const { return fTotEtAcc; }
+ Double_t GetTotEtAcc() const {
+ return fTotEtAcc;
+ }
- /** Total neutral Et in the event (without acceptance cuts) */
- Double_t GetTotNeutralEt() const { return fTotNeutralEt; }
+ /** Total neutral Et in the event (without acceptance cuts) */
+ Double_t GetTotNeutralEt() const {
+ return fTotNeutralEt;
+ }
/** Total neutral Et in the event within the acceptance cuts */
- Double_t GetTotNeutralEtAcc() const { return fTotNeutralEtAcc; }
-
+ Double_t GetTotNeutralEtAcc() const {
+ return fTotNeutralEtAcc;
+ }
+
/** Total charged Et in the event (without acceptance cuts) */
- Double_t GetTotChargedEt() const { return fTotChargedEt; }
+ Double_t GetTotChargedEt() const {
+ return fTotChargedEt;
+ }
/** Total charged Et in the event within the acceptance cuts */
- Double_t GetTotChargedEtAcc() const { return fTotChargedEtAcc; }
+ Double_t GetTotChargedEtAcc() const {
+ return fTotChargedEtAcc;
+ }
+
+ void SetTPCOnlyTrackCuts(const AliESDtrackCuts *cuts) {
+ fEsdtrackCutsTPC = (AliESDtrackCuts *) cuts;
+ }
- void SetTPCOnlyTrackCuts(const AliESDtrackCuts *cuts){ fEsdtrackCutsTPC = (AliESDtrackCuts *) cuts;}
-
/** Set the centrality object */
- void SetCentralityObject(AliCentrality *cent) { fCentrality = cent; }
-
+ void SetCentralityObject(AliCentrality *cent) {
+ fCentrality = cent;
+ }
+
/** Get contribution from non-removed charged particles */
- virtual Double_t GetChargedContribution(Int_t /*clusterMultiplicity*/) {return 0;}
+ virtual Double_t GetChargedContribution(Int_t /*clusterMultiplicity*/) {
+ return 0;
+ }
/** Get contribution from non-removed neutral particles */
- virtual Double_t GetNeutralContribution(Int_t /*clusterMultiplicity*/) {return 0;}
-
+ virtual Double_t GetNeutralContribution(Int_t /*clusterMultiplicity*/) {
+ return 0;
+ }
+
/** Get contribution from removed gammas */
- virtual Double_t GetGammaContribution(Int_t /*clusterMultiplicity*/) {return 0;}
+ virtual Double_t GetGammaContribution(Int_t /*clusterMultiplicity*/) {
+ return 0;
+ }
- void MakeSparseHistograms(){fMakeSparse=kTRUE;}
+ void MakeSparseHistograms() {
+ fMakeSparse=kTRUE;
+ }
+
+ AliAnalysisEtCuts * GetCuts() const { return fCuts; }
protected:
//AliAnalysisEtCuts *fCuts; // keeper of basic cuts
Double_t CalculateTransverseEnergy(AliESDCaloCluster *cluster);
- Double_t fTotEt;/** Total Et in the event (without acceptance cuts) */
- Double_t fTotEtAcc;/** Total Et in the event within the acceptance cuts */
+ Double_t TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, AliVEvent *e) const;
- Double_t fTotNeutralEt;/** Total neutral Et in the event */
- Double_t fTotNeutralEtAcc;/** Total neutral Et in the event within the acceptance cuts */
- Double_t fTotChargedEt;/** Total charged Et in the event */
+ virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &/*cluster*/) const { return kFALSE; }
+
+ Double_t fTotEt;/** Total Et in the event (without acceptance cuts) */
+ Double_t fTotEtAcc;/** Total Et in the event within the acceptance cuts */
+
+ Double_t fTotNeutralEt;/** Total neutral Et in the event */
+ Double_t fTotNeutralEtAcc;/** Total neutral Et in the event within the acceptance cuts */
+ Double_t fTotChargedEt;/** Total charged Et in the event */
Double_t fTotChargedEtAcc;/** Total charged Et in the event within the acceptance cuts */
- Int_t fMultiplicity;/** Multiplicity of particles in the event */
- Int_t fChargedMultiplicity;/** Multiplicity of charged particles in the event */
+ Int_t fMultiplicity;/** Multiplicity of particles in the event */
+ Int_t fChargedMultiplicity;/** Multiplicity of charged particles in the event */
Int_t fNeutralMultiplicity; /** Multiplicity of neutral particles in the event */
-
- Double_t fBaryonEt; /** Et of identified baryons; calo based (Rec only for now) */
+
+ Double_t fBaryonEt; /** Et of identified baryons; calo based (Rec only for now) */
Double_t fAntiBaryonEt; /** Et of identified anti-baryons; calo based (Rec only for now) */
Double_t fMesonEt; /** Et of identified mesons; calo based (Rec only for now) */
Double_t fNeutronEt; /** Et of neutrons (MC only for now) */
Double_t fAntiNeutronEt; /** Et of anti-neutrons (MC only for now) */
Double_t fGammaEt; /** Et of identified electrons (MC only for now) */
-
- Double_t fProtonEtAcc; /** Et of identified protons in calorimeter acceptance */
- Double_t fPionEtAcc; /** Et of identified pions in calorimeter acceptance */
- Double_t fChargedKaonEtAcc; /** Et of identified charged kaons in calorimeter acceptance */
+
+ Double_t fProtonEtAcc; /** Et of identified protons in calorimeter acceptance */
+ Double_t fPionEtAcc; /** Et of identified pions in calorimeter acceptance */
+ Double_t fChargedKaonEtAcc; /** Et of identified charged kaons in calorimeter acceptance */
Double_t fMuonEtAcc; /** Et of identified muons in calorimeter acceptance */
Double_t fElectronEtAcc; /** Et of identified electrons in calorimeter acceptance */
-
+
Float_t fEnergyDeposited; /** Energy deposited in calorimeter */
Float_t fEnergyTPC; /** Energy measured in TPC */
Short_t fCharge; /** Charge of the particle */
Bool_t fTrackPassedCut; /** The track is accepted by ESDTrackCuts */
Int_t fCentClass; // centrality class
-
+
Double_t fEtaCut;/** Cut in eta (standard |eta| < 0.5 )*/
/** Eta cut for our acceptance */
Double_t fEtaCutAcc; // Eta cut for our acceptance
-
- /** Min phi cut for our acceptance in radians */
- Double_t fPhiCutAccMin; // Min phi cut for our acceptance in radians
-
+
+ /** Min phi cut for our acceptance in radians */
+ Double_t fPhiCutAccMin; // Min phi cut for our acceptance in radians
+
/** Max phi cut for our acceptance in radians */
- Double_t fPhiCutAccMax; // Max phi cut for our acceptance in radians
-
+ Double_t fPhiCutAccMax; // Max phi cut for our acceptance in radians
+
/** Detector radius */
- Double_t fDetectorRadius; // Detector radius
-
- /** Cut on the cluster energy */
- Double_t fClusterEnergyCut; // Cut on the cluster energy
-
+ Double_t fDetectorRadius; // Detector radius
+
+ /** Cut on the cluster energy */
+ Double_t fClusterEnergyCut; // Cut on the cluster energy
+
/** Minimum energy to cut on single cell cluster */
Double_t fSingleCellEnergyCut; // Minimum energy to cut on single cell cluster
-
- Double_t fTrackDistanceCut; // cut on track distance
-
+
+ Double_t fTrackDistanceCut; // cut on track distance
+
Double_t fTrackDxCut; // cut on track distance in x
-
+
Double_t fTrackDzCut; // cut on track distance in z
+ Double_t fChargedEnergyRemoved;
+ Double_t fNeutralEnergyRemoved;
+ Double_t fGammaEnergyAdded;
+
// Declare the histograms
TH1F *fHistEt; //Et spectrum
/** The full charged Et spectrum measured */
- TH1F *fHistChargedEt; //Charged Et spectrum
+ TH1F *fHistChargedEt; //Charged Et spectrum
/** The full neutral Et spectrum measured */
TH1F *fHistNeutralEt; //Neutral Et spectrum
TH2F *fHistPhivsPtNeg; //phi vs pT plot for negative tracks
/* PID plots */
- TH1F *fHistBaryonEt; /** Et of identified baryons */
+ TH1F *fHistBaryonEt; /** Et of identified baryons */
TH1F *fHistAntiBaryonEt; /** Et of identified anti-baryons */
TH1F *fHistMesonEt; /** Et of identified mesons */
TH1F *fHistNeutronEt; /** Et of neutrons (MC only for now) */
TH1F *fHistAntiNeutronEt; /** Et of anti-neutrons (MC only for now) */
TH1F *fHistGammaEt; /** Et of gammas (MC only for now) */
-
- TH1F *fHistProtonEtAcc; /** Et of identified protons in calorimeter acceptance */
- TH1F *fHistPionEtAcc; /** Et of identified protons in calorimeter acceptance */
- TH1F *fHistChargedKaonEtAcc; /** Et of identified charged kaons in calorimeter acceptance */
+
+ TH1F *fHistProtonEtAcc; /** Et of identified protons in calorimeter acceptance */
+ TH1F *fHistPionEtAcc; /** Et of identified protons in calorimeter acceptance */
+ TH1F *fHistChargedKaonEtAcc; /** Et of identified charged kaons in calorimeter acceptance */
TH1F *fHistMuonEtAcc; /** Et of identified muons in calorimeter acceptance */
TH1F *fHistElectronEtAcc; /** Et of identified electrons in calorimeter acceptance */
-
+
/* Correction plots */
TH1F *fHistTMDeltaR; /* Track matching plots; Rec only for now */
TH2F *fHistTMDxDz; /* Track matching plots; Rec only for now */
-
+
/* Auxiliary Histogram variables */
static Float_t fgEtaAxis[17];//bins for eta axis of histograms
static Int_t fgnumOfEtaBins;//number of eta bins
static Int_t fgNumOfEBins;//number of pt bins
static Float_t fgRAxis[48];//bins for R axis
static Int_t fgNumOfRBins;//number of R bins
-
-
+
+
TTree *fTree; // optional TTree
TTree *fTreeDeposit; // optional TTree for energy deposit measurements
- /** Centrality object */
+ /** Centrality object */
AliCentrality *fCentrality; //Centrality object
-
+
Short_t fDetector; /** Which detector? (-1 -> PHOS, 1 -> EMCAL)*/
Bool_t fMakeSparse;//Boolean for whether or not to make sparse histograms
-
+
THnSparseF *fSparseHistTracks; /** THnSparse histograms */
-
+
THnSparseF *fSparseHistClusters; /** THnSparse histograms */
-
+
/** ET sparse valuses */
THnSparseF *fSparseHistEt; //!
-
+
/** Values for sparse hists */
Double_t *fSparseTracks; //!
-
+
/** Values for sparse hists */
Double_t *fSparseClusters; //!
-
+
/** ET sparse valuses */
Double_t *fSparseEt; //!
+
+ Char_t fClusterType; // selection on cluster type
+ Bool_t fMatrixInitialized;
+ TH1I *fCutFlow; // Cut flow
+ AliAnalysisEtSelector *fSelector;
private:
//Declare private to avoid compilation warning
TNamed()
//
,fCommonEtaCut(0.5)
- ,fCommonClusterEnergyCut(0.1)
+ ,fCommonClusterEnergyCut(0.15)
,fCommonTrackPtCut(0.0)
,fCommonSingleCell(1)
,fEmcalTrackDistanceCut(15.0)
,fPhosTrackDistanceCut(10.0)
,fPhosTrackDxCut(8.0)
,fPhosTrackDzCut(3.0)
+ ,fPhosTrackRCut(2.0)
+ ,fPhosBadDistanceCut(3.0)
,fGeometryPhosEtaAccCut(0.12)
,fGeometryPhosPhiAccMinCut(260.0)
,fReconstructedPidCut(0.0)
//
,fReconstructedPhosClusterType(-1)
- ,fReconstructedPhosClusterEnergyCut(0.1)
+ ,fReconstructedPhosClusterEnergyCut(0.15)
,fReconstructedPhosSingleCellEnergyCut(0.5)
,fReconstructedPhosTrackDistanceTightCut(3.0)
,fReconstructedPhosTrackDistanceMediumCut(5.0)
Int_t GetReconstructedNItsClustersCut() const { return fReconstructedNItsClustersCut; }
Double_t GetReconstructedPidCut() const { return fReconstructedPidCut; }
// ReconstructedPhos
- Char_t GetReconstructedPhosClusterType() const { return fReconstructedPhosClusterType; }
+ Char_t GetPhosClusterType() const { return fReconstructedPhosClusterType; }
Double_t GetReconstructedPhosClusterEnergyCut() const { return fReconstructedPhosClusterEnergyCut; }
Double_t GetReconstructedPhosSingleCellEnergyCut() const { return fReconstructedPhosSingleCellEnergyCut; }
Double_t GetPhosTrackDistanceCut() const { return fPhosTrackDistanceCut; }
Double_t GetPhosTrackDxCut() const { return fPhosTrackDxCut; }
Double_t GetPhosTrackDzCut() const { return fPhosTrackDzCut; }
+ Double_t GetPhosTrackRCut() const { return fPhosTrackRCut; }
+
+ Double_t GetPhosBadDistanceCut() const { return fPhosBadDistanceCut; }
+
// ReconstructedEmcal
- Char_t GetReconstructedEmcalClusterType() const { return fReconstructedEmcalClusterType; }
+ Char_t GetEmcalClusterType() const { return fReconstructedEmcalClusterType; }
Double_t GetReconstructedEmcalClusterEnergyCut() const { return fReconstructedEmcalClusterEnergyCut; }
Double_t GetReconstructedEmcalSingleCellEnergyCut() const { return fReconstructedEmcalSingleCellEnergyCut; }
Double_t GetEmcalTrackDistanceCut() const { return fEmcalTrackDistanceCut; }
Double_t GetHistMinParticlePt() const { return fHistMinParticlePt; }
Double_t GetHistMaxParticlePt() const { return fHistMaxParticlePt; }
+
+
Short_t GetDetectorPhos() const { return fgkDetectorPhos; }
Short_t GetDetectorEmcal() const { return fgkDetectorEmcal; }
void SetPhosTrackDistanceCut(const Double_t val) { fPhosTrackDistanceCut = val; }
void SetPhosTrackDxCut(const Double_t val) { fPhosTrackDxCut = val; }
void SetPhosTrackDzCut(const Double_t val) { fPhosTrackDzCut = val; }
+ void SetPhosTrackRCut(const Double_t val) { fPhosTrackRCut = val; }
+
+ void SetPhosBadDistanceCut(const Double_t val) { fPhosBadDistanceCut = val; }
// ReconstructedEmcal
void SetReconstructedEmcalClusterType(const Char_t val) { fReconstructedEmcalClusterType = val; }
Double_t fPhosTrackDistanceCut; // PHOS track distance
Double_t fPhosTrackDxCut; // PHOS track distance in x
Double_t fPhosTrackDzCut; // PHOS track distance in z
-
+ Double_t fPhosTrackRCut; // PHOS track distance in r (using the parametrized track distance)
+
+ Double_t fPhosBadDistanceCut; // PHOS distance to bad channel
+
// GeometryPhos
Double_t fGeometryPhosEtaAccCut; // PHOS Eta Acc cut
Double_t fGeometryPhosPhiAccMinCut; // PHOS Phi Acc Min cut
#include "AliLog.h"
#include <iostream>
#include <AliCentrality.h>
-
+#include "AliPHOSGeoUtils.h"
+#include "AliPHOSGeometry.h"
+#include "TFile.h"
using namespace std;
ClassImp(AliAnalysisEtMonteCarlo);
,fHistEtNonRemovedKaonPlus(0)
,fHistEtNonRemovedKaonMinus(0)
,fHistEtNonRemovedK0s(0)
+ ,fHistEtNonRemovedK0L(0)
,fHistEtNonRemovedLambdas(0)
,fHistEtNonRemovedElectrons(0)
,fHistEtNonRemovedPositrons(0)
,fHistEtNonRemovedAntiNeutrons(0)
,fHistEtNonRemovedGammas(0)
,fHistEtNonRemovedGammasFromPi0(0)
- ,fHistEtRemovedGammas(0)
+
+ ,fHistEtRemovedGammas(0)
,fHistEtRemovedNeutrons(0)
,fHistEtRemovedAntiNeutrons(0)
+
+ ,fHistEtRemovedCharged(0)
+ ,fHistEtRemovedNeutrals(0)
+ ,fHistEtNonRemovedCharged(0)
+ ,fHistEtNonRemovedNeutrals(0)
+
,fHistMultNonRemovedProtons(0)
,fHistMultNonRemovedAntiProtons(0)
,fHistMultNonRemovedPiPlus(0)
,fHistMultNonRemovedKaonPlus(0)
,fHistMultNonRemovedKaonMinus(0)
,fHistMultNonRemovedK0s(0)
+ ,fHistMultNonRemovedK0L(0)
,fHistMultNonRemovedLambdas(0)
,fHistMultNonRemovedElectrons(0)
,fHistMultNonRemovedPositrons(0)
,fHistMultNonRemovedNeutrons(0)
,fHistMultNonRemovedAntiNeutrons(0)
,fHistMultNonRemovedGammas(0)
+
,fHistMultRemovedGammas(0)
,fHistMultRemovedNeutrons(0)
,fHistMultRemovedAntiNeutrons(0)
+
+ ,fHistMultRemovedCharged(0)
+ ,fHistMultRemovedNeutrals(0)
+
+ ,fHistMultNonRemovedCharged(0)
+ ,fHistMultNonRemovedNeutrals(0)
+
,fHistTrackMultvsNonRemovedCharged(0)
,fHistTrackMultvsNonRemovedNeutral(0)
,fHistTrackMultvsRemovedGamma(0)
,fEtNonRemovedPiMinus(0)
,fEtNonRemovedKaonPlus(0)
,fEtNonRemovedKaonMinus(0)
- ,fEtNonRemovedK0s(0)
+ ,fEtNonRemovedK0S(0)
+ ,fEtNonRemovedK0L(0)
,fEtNonRemovedLambdas(0)
,fEtNonRemovedElectrons(0)
,fEtNonRemovedPositrons(0)
,fEtNonRemovedGammasFromPi0(0)
,fEtNonRemovedNeutrons(0)
,fEtNonRemovedAntiNeutrons(0)
+
+ ,fEtRemovedProtons(0)
+ ,fEtRemovedAntiProtons(0)
+ ,fEtRemovedPiPlus(0)
+ ,fEtRemovedPiMinus(0)
+ ,fEtRemovedKaonPlus(0)
+ ,fEtRemovedKaonMinus(0)
+ ,fEtRemovedK0s(0)
+ ,fEtRemovedK0L(0)
+ ,fEtRemovedLambdas(0)
+ ,fEtRemovedElectrons(0)
+ ,fEtRemovedPositrons(0)
+ ,fEtRemovedMuMinus(0)
+ ,fEtRemovedMuPlus(0)
+
+ ,fEtRemovedGammasFromPi0(0)
,fEtRemovedGammas(0)
,fEtRemovedNeutrons(0)
,fEtRemovedAntiNeutrons(0)
,fMultNonRemovedKaonPlus(0)
,fMultNonRemovedKaonMinus(0)
,fMultNonRemovedK0s(0)
+ ,fMultNonRemovedK0L(0)
,fMultNonRemovedLambdas(0)
,fMultNonRemovedElectrons(0)
,fMultNonRemovedPositrons(0)
,fMultNonRemovedGammas(0)
,fMultNonRemovedNeutrons(0)
,fMultNonRemovedAntiNeutrons(0)
+
+ ,fMultRemovedProtons(0)
+ ,fMultRemovedAntiProtons(0)
+ ,fMultRemovedPiPlus(0)
+ ,fMultRemovedPiMinus(0)
+ ,fMultRemovedKaonPlus(0)
+ ,fMultRemovedKaonMinus(0)
+ ,fMultRemovedK0s(0)
+ ,fMultRemovedK0L(0)
+ ,fMultRemovedLambdas(0)
+ ,fMultRemovedElectrons(0)
+ ,fMultRemovedPositrons(0)
+ ,fMultRemovedMuMinus(0)
+ ,fMultRemovedMuPlus(0)
+
+
,fMultRemovedGammas(0)
,fMultRemovedNeutrons(0)
,fMultRemovedAntiNeutrons(0)
+
,fTrackMultInAcc(0)
,fHistDxDzNonRemovedCharged(0)
,fHistDxDzRemovedCharged(0)
,fEnergyChargedRemoved(0)
,fEnergyChargedNotRemoved(0)
,fEnergyNeutralNotRemoved(0)
+ ,fNClusters(0)
+ ,fGeoUtils(0)
+ ,fBadMapM2(0)
+ ,fBadMapM3(0)
+ ,fBadMapM4(0)
{
fTrackDistanceCut = 7.0;
delete fHistEtNonRemovedKaonPlus; // enter comment here
delete fHistEtNonRemovedKaonMinus; // enter comment here
delete fHistEtNonRemovedK0s; // enter comment here
+ delete fHistEtNonRemovedK0L; // enter comment here
delete fHistEtNonRemovedLambdas; // enter comment here
delete fHistEtNonRemovedElectrons; // enter comment here
delete fHistEtNonRemovedPositrons; // enter comment here
delete fHistMultNonRemovedKaonPlus; // enter comment here
delete fHistMultNonRemovedKaonMinus; // enter comment here
delete fHistMultNonRemovedK0s; // enter comment here
+ delete fHistMultNonRemovedK0L; // enter comment here
delete fHistMultNonRemovedLambdas; // enter comment here
delete fHistMultNonRemovedElectrons; // enter comment here
delete fHistMultNonRemovedPositrons; // enter comment here
fImpactParameter = hijingGenHeader->ImpactParameter();
fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
- /*
- printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
- hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
- printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
- hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants());
- printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
- hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
- printf("NN %d NNw %d NwN %d, NwNw %d\n",
- hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
- */
- }
+ }
/* // placeholder if we want to get some Pythia info later
AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
}
*/
-
// Let's play with the stack!
AliStack *stack = event->Stack();
// Int_t iPartLastMom = part->GetMother(1);
TParticlePDG *pdg = part->GetPDG(0);
- //TParticlePDG *pdgMom = 0;
-// TParticlePDG *pdgMomLast = 0;
- if (!pdg)
+ if (!pdg)
{
//Printf("ERROR: Could not get particle PDG %d", iPart);
continue;
}
-// if (iPartMom>0)
-// {
-// partMom = stack->Particle(iPartMom);
-// //pdgMom = partMom->GetPDG(0);
-// }
-
-// if (iPartLastMom>0)
-// {
-// partMomLast = stack->Particle(iPartLastMom);
-// pdgMomLast = partMomLast->GetPDG(0);
-// }
-
-
Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
-
- // Check if it is a primary particle
- //if (!stack->IsPhysicalPrimary(iPart)) continue;
-
- //printf("MC: iPart %03d eta %4.3f phi %4.3f code %d charge %g \n", iPart, part->Eta(), part->Phi(), pdg->PdgCode(), pdg->Charge()); // tmp/debug printout
+ Int_t code = pdg->PdgCode();
// Check for reasonable (for now neutral and singly charged) charge on the particle
//TODO:Maybe not only singly charged?
{
// calculate E_T
if (
- TMath::Abs(pdg->PdgCode()) == fgProtonCode ||
- TMath::Abs(pdg->PdgCode()) == fgNeutronCode ||
- TMath::Abs(pdg->PdgCode()) == fgLambdaCode ||
- TMath::Abs(pdg->PdgCode()) == fgXiCode ||
- TMath::Abs(pdg->PdgCode()) == fgXi0Code ||
- TMath::Abs(pdg->PdgCode()) == fgOmegaCode
+ TMath::Abs(code) == fgProtonCode ||
+ TMath::Abs(code) == fgNeutronCode ||
+ TMath::Abs(code) == fgLambdaCode ||
+ TMath::Abs(code) == fgXiCode ||
+ TMath::Abs(code) == fgXi0Code ||
+ TMath::Abs(code) == fgOmegaCode
)
{
- if (pdg->PdgCode() > 0) {
+ if (code > 0) {
particleMassPart = - protonMass;
}
- if (pdg->PdgCode() < 0) {
+ if (code < 0) {
particleMassPart = protonMass;
}
}
Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
- fSparseTracks[0] = pdg->PdgCode();
+ fSparseTracks[0] = code;
fSparseTracks[1] = pdg->Charge()*3;
fSparseTracks[2] = pdg->Mass();
fSparseTracks[3] = et;
fSparseHistTracks->Fill(fSparseTracks);
// Fill up total E_T counters for each particle species
- if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
+ if (code == fgProtonCode || code == fgAntiProtonCode)
{
fProtonEt += et;
}
- if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
+ if (code == fgPiPlusCode || code == fgPiMinusCode)
{
fPionEt += et;
- if (pdg->PdgCode() == fgPiPlusCode)
+ if (code == fgPiPlusCode)
{
fPiPlusMult++;
}
fPiMinusMult++;
}
}
- if (pdg->PdgCode() == fgGammaCode)
+ if (code == fgGammaCode)
{
fPiZeroMult++;
}
- if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
+ if (code == fgKPlusCode || code == fgKMinusCode)
{
fChargedKaonEt += et;
}
- if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
+ if (code == fgMuPlusCode || code == fgMuMinusCode)
{
fMuonEt += et;
}
- if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+ if (code == fgEPlusCode || code == fgEMinusCode)
{
fElectronEt += et;
}
-
// some neutrals also
- if (pdg->PdgCode() == fgNeutronCode)
+ if (code == fgNeutronCode)
{
fNeutronEt += et;
}
- if (pdg->PdgCode() == fgAntiNeutronCode)
+ if (code == fgAntiNeutronCode)
{
fAntiNeutronEt += et;
}
- if (pdg->PdgCode() == fgGammaCode)
+ if (code == fgGammaCode)
{
fGammaEt += et;
}
// Neutral particles
- if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
+ //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
+
+ if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
{
if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
{
fNeutralMultiplicity++;
- //std::cout << pdg->PdgCode() << ", " << iPart << ", " << part->GetMother(0) << ", " << stack->IsPhysicalPrimary(iPart) << ", " << part->GetNDaughters() << ", " << part->Energy() << ", " << part->Eta() << ", " << part->Phi() << std::endl;
+ //std::cout << code << ", " << iPart << ", " << part->GetMother(0) << ", " << stack->IsPhysicalPrimary(iPart) << ", " << part->GetNDaughters() << ", " << part->Energy() << ", " << part->Eta() << ", " << part->Phi() << std::endl;
//if(part->GetDaughter(0) > 0) std::cout << stack->Particle(part->GetDaughter(0))->GetPdgCode() << " " << stack->Particle(part->GetDaughter(1))->GetPdgCode() << std::endl;
- if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEtAcc += et;
- if (et > fCuts->GetCommonClusterEnergyCut()) fTotEtAcc += et;
+ fTotNeutralEtAcc += et;
+ fTotEtAcc += et;
+ //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEtAcc += et;
+ //if (et > fCuts->GetCommonClusterEnergyCut()) fTotEtAcc += et;
+
if (part->Energy() > 0.05) partCount++;
}
}
fTotChargedEtAcc += et;
fTotEtAcc += et;
- if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
+ if (code == fgProtonCode || code == fgAntiProtonCode)
{
fProtonEtAcc += et;
}
- if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
+ if (code == fgPiPlusCode || code == fgPiMinusCode)
{
fPionEtAcc += et;
}
- if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
+ if (code == fgKPlusCode || code == fgKMinusCode)
{
fChargedKaonEtAcc += et;
}
- if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
+ if (code == fgMuPlusCode || code == fgMuMinusCode)
{
fMuonEtAcc += et;
}
- if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+ if (code == fgEPlusCode || code == fgEMinusCode)
{
fElectronEtAcc += et;
}
}
fTotEt = fTotChargedEt + fTotNeutralEt;
- fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
-
-
+ fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
//std::cout << "Event done! # of particles: " << partCount << std::endl;
return 0;
}
AliFatal("ERROR: Event does not exist");
return 0;
}
+ //if (!fMatrixInitialized)
+ if(0)
+ {
+ for (Int_t mod=0; mod<5; mod++) {
+ if (!ev->GetPHOSMatrix(mod)) continue;
+ fGeoUtils->SetMisalMatrix(ev->GetPHOSMatrix(mod),mod) ;
+ Printf("PHOS geo matrix %p for module # %d is set\n", ev->GetPHOSMatrix(mod), mod);
+ }
+ fMatrixInitialized = kTRUE;
+ }
fSparseClusters[0] = 0;
fSparseClusters[1] = 0;
fSparseClusters[2] = 0;
AliFatal("Detector ID has not been specified");
return -1;
}
-
+if(0)
+{
// Track loop to fill a pT spectrum
for (Int_t iTracks = 0; iTracks < realEvent->GetNumberOfTracks(); iTracks++)
{
fTrackMultInAcc++;
}
}
-
+}
Int_t nCluster = caloClusters->GetEntries();
-
+ fNClusters = 0;
// loop the clusters
for (int iCluster = 0; iCluster < nCluster; iCluster++ )
{
+
AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
//Float_t caloE = caloCluster->E();
+
+ if(caloCluster->GetType() != fClusterType) continue;
+ if(AliAnalysisEtMonteCarlo::TooCloseToBadChannel(*caloCluster))continue;
+ fNClusters++;
UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
- TParticle *part = stack->Particle(iPart);
- TParticle *partMom = 0;
+ TParticle *part = stack->Particle(iPart);
+ //TParticle *partMom = 0;
if (!part)
{
// compare MC and Rec energies for all particles
//fHistAllERecEMC->Fill(part->Energy(),caloE);
- TParticlePDG *pdg = part->GetPDG(0);
- TParticlePDG *pdgMom = 0;
-
- if (!pdg)
- {
- Printf("ERROR: Could not get particle PDG %d", iPart);
- continue;
- }
+ //TParticlePDG *pdgMom = 0;
+/*
Int_t iPartMom = part->GetMother(0);
if (iPartMom>0)
{
partMom = stack->Particle(iPartMom);
pdgMom = partMom->GetPDG(0);
- }
-
- // Check if it is a primary particle
- //if (!stack->IsPhysicalPrimary(iPart)) continue;
+ }*/
+ //int mcClusterCharge = stack->Particle(iPart)->GetPDG()->Charge();
if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
{
- if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
- continue;
+ // if (!((code == fgEPlusCode) || (code == fgEMinusCode) || (code == fgGammaCode)))
+ int p = GetPrimMother(iPart, stack);
+ // std::cout << p << std::endl;
+ if(p > 0)
+ {
+ part = stack->Particle(p);
+ }
+ else
+ {
+ std::cout << "What!? No primary?" << std::endl;
+ continue;
+ }
+
} // end of primary particle check
-
- if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
+ TParticlePDG *pdg = part->GetPDG();
+ Int_t code = pdg->PdgCode();
+ if (!pdg)
+ {
+ Printf("ERROR: Could not get particle PDG %d", iPart);
+ continue;
+ }
+ //if(code == fgK0SCode)PrintFamilyTree(iPart, stack);
+ //pdgMom = stack->Particle(GetPrimMother(iPart, stack))->GetPDG();
+
+ //-if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
Double_t clEt = CalculateTransverseEnergy(caloCluster);
+// if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
if (caloCluster->E() < fCuts->GetCommonClusterEnergyCut()) continue;
Float_t pos[3];
caloCluster->GetPosition(pos);
TVector3 cp(pos);
- fSparseClusters[0] = pdg->PdgCode();
+ fSparseClusters[0] = code;
fSparseClusters[1] = pdg->Charge()/3;
fSparseClusters[2] = pdg->Mass();
fSparseClusters[3] = clEt;
fSparseClusters[10] = 0;
fSparseHistClusters->Fill(fSparseClusters);
+
//if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
-
- if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
+ pdg = part->GetPDG(0);
+ //if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
+ if(caloCluster->GetTrackMatchedIndex() > -1 && (fCuts->GetPhosTrackRCut() > TestCPV(caloCluster->GetTrackDx(), caloCluster->GetTrackDz(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Pt(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Charge(), ev)))
{
if (pdg->Charge() != 0)
{
- //std::cout << "Removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+ //std::cout << "Removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
fChargedRemoved++;
- fEnergyChargedRemoved += caloCluster->E();
+ //fEnergyChargedRemoved += caloCluster->E();
+ fEnergyChargedRemoved += clEt;
fHistRemovedOrNot->Fill(0.0, fCentClass);
//std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
//fHistDecayVertexRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
}
else
{
- //std::cout << "Removing neutral: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+ //std::cout << "Removing neutral: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
//std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
fNeutralRemoved++;
- fEnergyNeutralRemoved += caloCluster->E();
+ //fEnergyNeutralRemoved += caloCluster->E();
+ fEnergyNeutralRemoved += clEt;
fHistRemovedOrNot->Fill(1.0, fCentClass);
//std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
//fHistDecayVertexRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
- if (pdg->PdgCode() == fgGammaCode)
+ if (code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
{
fEtRemovedGammas += clEt;
fMultRemovedGammas++;
}
- else if (pdg->PdgCode() == fgNeutronCode)
+ else if (code == fgNeutronCode)
{
fEtRemovedNeutrons += clEt;
fMultRemovedNeutrons++;
}
- else if (pdg->PdgCode() == fgAntiNeutronCode)
+ else if (code == fgAntiNeutronCode)
{
fEtRemovedAntiNeutrons += clEt;
fMultRemovedAntiNeutrons++;
}
- //else std::cout << "Hmm, what is this (neutral: " << pdg->PdgCode() << std::endl;
+
+ //else std::cout << "Hmm, what is this (neutral: " << code << std::endl;
}
}
else
if (pdg->Charge() != 0)
{
- //std::cout << "Not removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+ std::cout << "Not removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
//std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
fChargedNotRemoved++;
- fEnergyChargedNotRemoved += caloCluster->E();
+ //fEnergyChargedNotRemoved += caloCluster->E();
+ fEnergyChargedNotRemoved += clEt;
fHistRemovedOrNot->Fill(2.0, fCentClass);
//std::cout << fHistDecayVertexNonRemovedCharged << std::endl;
//std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
//fHistDecayVertexNonRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
- if (pdg->PdgCode() == fgProtonCode)
+ if (code == fgProtonCode)
{
//std::cout << clEt << std::endl;
fEtNonRemovedProtons += clEt;
fMultNonRemovedProtons++;
}
- else if (pdg->PdgCode() == fgAntiProtonCode)
+ else if (code == fgAntiProtonCode)
{
//std::cout << clEt << std::endl;
fEtNonRemovedAntiProtons += clEt;
fMultNonRemovedAntiProtons++;
}
- else if (pdg->PdgCode() == fgPiPlusCode)
+ else if (code == fgPiPlusCode)
{
//std::cout << "PI+" << clEt << std::endl;
fEtNonRemovedPiPlus += clEt;
fMultNonRemovedPiPlus++;
}
- else if (pdg->PdgCode() == fgPiMinusCode)
+ else if (code == fgPiMinusCode)
{
// std::cout << "PI-" << clEt << std::endl;
fEtNonRemovedPiMinus += clEt;
fMultNonRemovedPiMinus++;
}
- else if (pdg->PdgCode() == fgKPlusCode)
+ else if (code == fgKPlusCode)
{
//std::cout << clEt << std::endl;
fEtNonRemovedKaonPlus += clEt;
fMultNonRemovedKaonPlus++;
}
- else if (pdg->PdgCode() == fgKMinusCode)
+ else if (code == fgKMinusCode)
{
//std::cout << clEt << std::endl;
fEtNonRemovedKaonMinus += clEt;
fMultNonRemovedKaonMinus++;
}
- else if (pdg->PdgCode() == fgEPlusCode)
+ else if (code == fgEPlusCode)
{
//std::cout << clEt << std::endl;
if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
fMultNonRemovedPositrons++;
}
}
- else if (pdg->PdgCode() == fgEMinusCode)
+ else if (code == fgEMinusCode)
{
//std::cout << clEt << std::endl;
if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
fMultNonRemovedElectrons++;
}
}
- else if (pdg->PdgCode() == fgMuPlusCode)
+ else if (code == fgMuPlusCode)
{
- std::cout << clEt << std::endl;
fEtNonRemovedMuPlus += clEt;
fMultNonRemovedMuPlus++;
}
- else if (pdg->PdgCode() == fgMuMinusCode)
+ else if (code == fgMuMinusCode)
{
- std::cout << clEt << std::endl;
fEtNonRemovedMuMinus += clEt;
fMultNonRemovedMuMinus++;
}
}
else
{
- //std::cout << "Accepted: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
- //std::cout << "Not removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+ //std::cout << "Accepted: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+ //std::cout << "Not removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
//std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << stack->Particle(part->GetMother(0))->GetPDG()->GetName() << ", E: " << part->Energy() << std::endl;
fNeutralNotRemoved++;
- fEnergyNeutralNotRemoved += caloCluster->E();
+ fEnergyNeutralNotRemoved += clEt;
fHistRemovedOrNot->Fill(3.0, fCentClass);
- //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
- //fHistDecayVertexNonRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
- if (pdg->PdgCode() == fgGammaCode)
+ if (code == fgGammaCode)
{
fEtNonRemovedGammas += clEt;
fMultNonRemovedGammas++;
- if (pdgMom)
- {
- if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
- {
+// if (pdgMom)
+ // {
+ // if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
+ // {
// std::cout << "Mother of gamma: " << pdgMom->PdgCode() << " " << pdgMom->GetName() << ", E: " << part->Energy() << std::endl;
- fEtNonRemovedGammasFromPi0 += clEt;
- }
- }
+ // fEtNonRemovedGammasFromPi0 += clEt;
+ // }
+ // }
}
- else if (pdg->PdgCode() == fgNeutronCode)
+ else if(TMath::Abs(code) == fgPi0Code)
+ {
+ fEtNonRemovedGammasFromPi0 += clEt;
+ fMultNonRemovedGammas++;
+ }
+ else if(TMath::Abs(code) == fgEtaCode)
+ {
+ fEtNonRemovedGammasFromPi0 += clEt;
+ fMultNonRemovedGammas++;
+ }
+ else if(TMath::Abs(code) == 331)
+ {
+ fEtNonRemovedGammasFromPi0 += clEt;
+ fMultNonRemovedGammas++;
+ }
+ else if (code == fgNeutronCode)
{
fEtNonRemovedNeutrons += clEt;
fMultNonRemovedNeutrons++;
}
- else if (pdg->PdgCode() == fgAntiNeutronCode)
+ else if (code == fgAntiNeutronCode)
{
fEtNonRemovedAntiNeutrons += clEt;
fMultNonRemovedAntiNeutrons++;
}
- else if (pdg->PdgCode() == fgK0LCode || pdg->PdgCode() == fgK0SCode)
+ //else if (code == fgK0LCode || pdgMom->PdgCode() == fgK0SCode)
+ else if (code == fgK0SCode)
{
- fEtNonRemovedK0s += clEt;
+ //std::cout << "K0 with energy: " << clEt << std::endl;
+ fEtNonRemovedK0S += clEt;
fMultNonRemovedK0s++;
}
- else if (TMath::Abs(pdg->PdgCode()) == fgLambdaCode)
+ else if(TMath::Abs(code) == fgK0LCode)
+ {
+
+ fEtNonRemovedK0L += clEt;
+ fMultNonRemovedK0L++;
+ }
+
+ else if (TMath::Abs(code) == fgLambdaCode)
{
fEtNonRemovedLambdas += clEt;
fMultNonRemovedLambdas++;
}
- else std::cout << "Hmm, what is this (neutral not removed): " << pdg->PdgCode() << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
+
+ else std::cout << "Hmm, what is this (neutral not removed): " << code << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
}
}
} // end of loop over clusters
void AliAnalysisEtMonteCarlo::Init()
{ // init
AliAnalysisEt::Init();
+ TFile *f = TFile::Open("badchannels.root", "READ");
+
+ fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
+ fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
+ fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
+
}
void AliAnalysisEtMonteCarlo::ResetEventValues()
fEtNonRemovedPiMinus = 0;
fEtNonRemovedKaonPlus = 0;
fEtNonRemovedKaonMinus = 0;
- fEtNonRemovedK0s = 0;
+ fEtNonRemovedK0S = 0;
+ fEtNonRemovedK0L = 0;
fEtNonRemovedLambdas = 0;
fEtNonRemovedElectrons = 0;
fEtNonRemovedPositrons = 0;
fEtNonRemovedGammas = 0;
fEtNonRemovedGammasFromPi0 = 0;
+ fEtRemovedProtons = 0;
+ fEtRemovedAntiProtons = 0;
+ fEtRemovedPiPlus = 0;
+ fEtRemovedPiMinus = 0;
+ fEtRemovedKaonPlus = 0;
+ fEtRemovedKaonMinus = 0;
+ fEtRemovedK0s = 0;
+ fEtRemovedK0L = 0;
+ fEtRemovedLambdas = 0;
+ fEtRemovedElectrons = 0;
+ fEtRemovedPositrons = 0;
+ fEtRemovedMuPlus = 0;
+ fEtRemovedMuMinus = 0;
+ fEtRemovedNeutrons = 0;
+
+ fEtRemovedGammasFromPi0 = 0;
fEtRemovedGammas = 0;
fEtRemovedNeutrons = 0;
fEtRemovedAntiNeutrons = 0;
fMultNonRemovedKaonPlus = 0;
fMultNonRemovedKaonMinus = 0;
fMultNonRemovedK0s = 0;
+ fMultNonRemovedK0L = 0;
fMultNonRemovedLambdas = 0;
fMultNonRemovedElectrons = 0;
fMultNonRemovedPositrons = 0;
fMultNonRemovedAntiNeutrons = 0;
fMultNonRemovedGammas = 0;
+ fMultRemovedProtons = 0;
+ fMultRemovedAntiProtons = 0;
+ fMultRemovedPiPlus = 0;
+ fMultRemovedPiMinus = 0;
+ fMultRemovedKaonPlus = 0;
+ fMultRemovedKaonMinus = 0;
+ fMultRemovedK0s = 0;
+ fMultRemovedK0L = 0;
+ fMultRemovedLambdas = 0;
+ fMultRemovedElectrons = 0;
+ fMultRemovedPositrons = 0;
+ fMultRemovedMuPlus = 0;
+ fMultRemovedMuMinus = 0;
+
fMultRemovedGammas = 0;
fMultRemovedNeutrons = 0;
fMultRemovedAntiNeutrons = 0;
+ fEnergyChargedNotRemoved = 0;
+ fEnergyChargedRemoved = 0;
+ fEnergyNeutralNotRemoved = 0;
+ fEnergyNeutralRemoved = 0;
+
+ fChargedNotRemoved = 0;
+ fChargedRemoved = 0;
+ fNeutralNotRemoved = 0;
+ fNeutralRemoved = 0;
+
+
fTrackMultInAcc = 0;
}
fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
+ fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtRemovedGammas = new TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtRemovedNeutrons = new TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
fHistEtRemovedAntiNeutrons = new TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
+
+ fHistEtRemovedCharged = new TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
+ fHistEtRemovedNeutrals = new TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
+ fHistEtNonRemovedCharged = new TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
+ fHistEtNonRemovedNeutrals = new TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
+
fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
+ fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
- fHistMultNonRemovedGammas = new TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
+ fHistMultNonRemovedGammas = new TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
- fHistMultRemovedGammas = new TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
+ fHistMultRemovedGammas = new TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
fHistMultRemovedNeutrons = new TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
fHistMultRemovedAntiNeutrons = new TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
+ /*
+ fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
+ fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
+
+ fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
+ fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
+
+
+ fHistMultRemovedCharged = new TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
+ fHistMultRemovedNeutrals = new TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
+
+ fHistMultNonRemovedCharged = new TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
+ fHistMultNonRemovedNeutrals = new TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
list->Add(fHistEtNonRemovedKaonPlus);
list->Add(fHistEtNonRemovedKaonMinus);
list->Add(fHistEtNonRemovedK0s);
+ list->Add(fHistEtNonRemovedK0L);
list->Add(fHistEtNonRemovedLambdas);
list->Add(fHistEtNonRemovedElectrons);
list->Add(fHistEtNonRemovedPositrons);
list->Add(fHistEtRemovedNeutrons);
list->Add(fHistEtRemovedAntiNeutrons);
+ list->Add(fHistEtRemovedCharged);
+ list->Add(fHistEtRemovedNeutrals);
+
+ list->Add(fHistEtNonRemovedCharged);
+ list->Add(fHistEtNonRemovedNeutrals);
list->Add(fHistMultNonRemovedProtons);
list->Add(fHistMultNonRemovedAntiProtons);
list->Add(fHistMultNonRemovedKaonPlus);
list->Add(fHistMultNonRemovedKaonMinus);
list->Add(fHistMultNonRemovedK0s);
+ list->Add(fHistMultNonRemovedK0L);
list->Add(fHistMultNonRemovedLambdas);
list->Add(fHistMultNonRemovedElectrons);
list->Add(fHistMultNonRemovedPositrons);
list->Add(fHistMultRemovedNeutrons);
list->Add(fHistMultRemovedAntiNeutrons);
+ list->Add(fHistMultRemovedCharged);
+ list->Add(fHistMultRemovedNeutrals);
+
+ list->Add(fHistMultNonRemovedCharged);
+ list->Add(fHistMultNonRemovedNeutrals);
+
list->Add(fHistTrackMultvsNonRemovedCharged);
list->Add(fHistTrackMultvsNonRemovedNeutral);
list->Add(fHistTrackMultvsRemovedGamma);
fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
- fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0s, fCentClass);
+ fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
+ fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
- fHistEtRemovedGammas->Fill(fEtRemovedGammas, fCentClass);
+ fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
+// fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
+// +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
+// +fEtRemovedProtons.
+// fCentClass);
+// fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
+//
+// fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
+// +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
+// +fEtNonRemovedProtons,
+// fCentClass);
+// fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
+
+ fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
+ fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
+ fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
+ fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
+
+ fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
+ fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
+ fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
+ fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
+
+
fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
+ fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
+fMultNonRemovedProtons);
fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
- fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
+ fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
fMultRemovedGammas);
- fHistClusterMultvsNonRemovedCharged->Fill(fNeutralMultiplicity,
+ fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
+fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
- fHistClusterMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
- fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
+ fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
+ fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
- fHistClusterMultvsRemovedGamma->Fill(fTrackMultInAcc,
+ fHistClusterMultvsRemovedGamma->Fill(fNClusters,
fMultRemovedGammas);
}
+Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
+{
+ TParticle *part = stack->Particle(partIdx);
+ std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
+ std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
+ std::cout << "Energy: " << part->Energy() << std::endl;
+ std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
+ PrintMothers(partIdx, stack, 1);
+ return 0;
+}
+
+Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
+{
+ char *tabs = new char[gen+1];
+ for(Int_t i = 0; i < gen; ++i)
+ {
+ //std::cout << i << std::endl;
+ tabs[i] = '\t';
+ }
+ tabs[gen] = '\0';
+ Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
+ if(mothIdx < 0) return 0;
+ TParticle *mother = stack->Particle(mothIdx);
+ std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
+ std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
+ std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
+ std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
+ if(mother->GetPdgCode() == fgK0SCode)
+ {
+// std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
+ }
+// std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
+// std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
+// std::cout << "Energy: " << mother->Energy() << std::endl;
+// std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
+
+ delete [] tabs;
+ return PrintMothers(mothIdx, stack, gen+1) + 1;
+}
+
+Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
+{
+ if(partIdx >= 0)
+ {
+ Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
+ if(mothIdx < 0) return -1;
+ TParticle *mother = stack->Particle(mothIdx);
+ if(mother)
+ {
+ // if(mother->GetPdgCode() == fgK0SCode)
+ //{
+// std::cout << "!!!!!!!!!!!!!!!!! K0S !!!!!!!!!!!!!!!!!!" << std::endl;
+ //return mothIdx;
+ // }
+ //if(mother->GetPdgCode() == fgK0SCode&& stack->IsPhysicalPrimary(mothIdx))
+ //{
+// std::cout << "!!!!!!!!!!!!!!!!! Primary K0S !!!!!!!!!!!!!!!!!!" << std::endl;
+ //return mothIdx;
+ // }
+ if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
+ else return GetPrimMother(mothIdx, stack);
+ }
+ else
+ {
+ std::cout << "WAT!" << std::endl;
+ return -1;
+ }
+ }
+ return -1;
+}
+
+Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
+{
+ if(partIdx >= 0)
+ {
+ if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
+ Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
+ if(mothIdx < 0) return -1;
+ TParticle *mother = stack->Particle(mothIdx);
+ if(mother)
+ {
+ if(mother->GetPdgCode() == fgK0SCode)
+ {
+ return mothIdx;
+ }
+ return GetK0InFamily(mothIdx, stack);
+ }
+ else
+ {
+ std::cout << "WAT WAT!" << std::endl;
+ return -1;
+ }
+ }
+ return -1;
+}
+
+
+
+
+Bool_t AliAnalysisEtMonteCarlo::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
+{
+ Float_t gPos[3];
+ cluster.GetPosition(gPos);
+ Int_t relId[4];
+ TVector3 glVec(gPos);
+ fGeoUtils->GlobalPos2RelId(glVec, relId);
+
+ TVector3 locVec;
+ fGeoUtils->Global2Local(locVec, glVec, relId[0]);
+// std::cout << fGeoUtils << std::endl;
+ //std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
+ //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
+ for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
+ {
+ for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
+ {
+ if (relId[0] == 3)
+ {
+ if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 3;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 2)
+ {
+ if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 2;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 1)
+ {
+ if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 1;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kTRUE;
+ }
+ }
+ }
+
+ }
+ }
+ return kFALSE;
+
+
+}
#ifndef ALIANALYSISETMONTECARLO_H
#define ALIANALYSISETMONTECARLO_H
-//_________________________________________________________________________
+
+class AliPHOSGeometry;//_________________________________________________________________________
// Utility Class for transverse energy studies
// Base class for MC analysis
// - MC output
#include "AliAnalysisEt.h"
class TParticle;
class TH3F;
+class TH2I;
+class AliPHOSGeoUtils;
+class AliStack;
//class AliMCEvent;
//class AliESDEvent;
{
public:
-
- AliAnalysisEtMonteCarlo();
- virtual ~AliAnalysisEtMonteCarlo();
+
+ AliAnalysisEtMonteCarlo();
+ virtual ~AliAnalysisEtMonteCarlo();
virtual Int_t AnalyseEvent(AliVEvent* event);
- virtual Int_t AnalyseEvent(AliVEvent* event, AliVEvent* event2);
- //virtual Int_t AnalyseEvent(AliMCEvent* event, AliESDEvent* event2);
+ virtual Int_t AnalyseEvent(AliVEvent* event, AliVEvent* event2);
+ //virtual Int_t AnalyseEvent(AliMCEvent* event, AliESDEvent* event2);
virtual void Init();
virtual void ResetEventValues();
protected:
virtual bool TrackHitsCalorimeter(TParticle *part, Double_t magField=0.5);
+
+
+ Int_t GetPrimMother(Int_t partIdx, AliStack *stack);
+
+ Int_t GetK0InFamily(Int_t partIdx, AliStack *stack);
+
+ Int_t PrintFamilyTree(Int_t partIdx, AliStack *stack);
+ Int_t PrintMothers(Int_t partIdx, AliStack *stack, Int_t gen);
+
+
+ virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
+
protected:
Int_t fNcoll; // Ncoll, for Hijing; 1 otherwise
Int_t fNpart; // Ncoll, for Hijing; 2 otherwise
- TH3F *fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
- TH3F *fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
- TH3F *fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
- TH3F *fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
-
- TH2F *fHistRemovedOrNot; // If charged/neutral particles were removed or not
-
- TH2F *fHistEtNonRemovedProtons; // enter comment here
- TH2F *fHistEtNonRemovedAntiProtons; // enter comment here
- TH2F *fHistEtNonRemovedPiPlus; // enter comment here
- TH2F *fHistEtNonRemovedPiMinus; // enter comment here
- TH2F *fHistEtNonRemovedKaonPlus; // enter comment here
- TH2F *fHistEtNonRemovedKaonMinus; // enter comment here
- TH2F *fHistEtNonRemovedK0s; // enter comment here
- TH2F *fHistEtNonRemovedLambdas; // enter comment here
- TH2F *fHistEtNonRemovedElectrons; // enter comment here
- TH2F *fHistEtNonRemovedPositrons; // enter comment here
- TH2F *fHistEtNonRemovedMuPlus; // enter comment here
- TH2F *fHistEtNonRemovedMuMinus; // enter comment here
- TH2F *fHistEtNonRemovedNeutrons; // enter comment here
- TH2F *fHistEtNonRemovedAntiNeutrons; // enter comment here
- TH2F *fHistEtNonRemovedGammas; // enter comment here
- TH2F *fHistEtNonRemovedGammasFromPi0; // enter comment here
-
- TH2F *fHistEtRemovedGammas; // enter comment here
- TH2F *fHistEtRemovedNeutrons; // enter comment here
- TH2F *fHistEtRemovedAntiNeutrons; // enter comment here
-
-
- TH2F *fHistMultNonRemovedProtons; // enter comment here
- TH2F *fHistMultNonRemovedAntiProtons; // enter comment here
- TH2F *fHistMultNonRemovedPiPlus; // enter comment here
- TH2F *fHistMultNonRemovedPiMinus; // enter comment here
- TH2F *fHistMultNonRemovedKaonPlus; // enter comment here
- TH2F *fHistMultNonRemovedKaonMinus; // enter comment here
- TH2F *fHistMultNonRemovedK0s; // enter comment here
- TH2F *fHistMultNonRemovedLambdas; // enter comment here
- TH2F *fHistMultNonRemovedElectrons; // enter comment here
- TH2F *fHistMultNonRemovedPositrons; // enter comment here
- TH2F *fHistMultNonRemovedMuPlus; // enter comment here
- TH2F *fHistMultNonRemovedMuMinus; // enter comment here
- TH2F *fHistMultNonRemovedNeutrons; // enter comment here
- TH2F *fHistMultNonRemovedAntiNeutrons; // enter comment here
- TH2F *fHistMultNonRemovedGammas; // enter comment here
-
- TH2F *fHistMultRemovedGammas; // enter comment here
- TH2F *fHistMultRemovedNeutrons; // enter comment here
- TH2F *fHistMultRemovedAntiNeutrons; // enter comment here
-
- TH2F *fHistTrackMultvsNonRemovedCharged; // enter comment here
- TH2F *fHistTrackMultvsNonRemovedNeutral; // enter comment here
- TH2F *fHistTrackMultvsRemovedGamma; // enter comment here
-
- TH2F *fHistClusterMultvsNonRemovedCharged; // enter comment here
- TH2F *fHistClusterMultvsNonRemovedNeutral; // enter comment here
- TH2F *fHistClusterMultvsRemovedGamma; // enter comment here
-
- TH2F *fHistMultvsNonRemovedChargedE; // enter comment here
- TH2F *fHistMultvsNonRemovedNeutralE; // enter comment here
- TH2F *fHistMultvsRemovedGammaE; // enter comment here
-
- Float_t fEtNonRemovedProtons; // enter comment here
- Float_t fEtNonRemovedAntiProtons; // enter comment here
- Float_t fEtNonRemovedPiPlus; // enter comment here
- Float_t fEtNonRemovedPiMinus; // enter comment here
- Float_t fEtNonRemovedKaonPlus; // enter comment here
- Float_t fEtNonRemovedKaonMinus; // enter comment here
- Float_t fEtNonRemovedK0s; // enter comment here
- Float_t fEtNonRemovedLambdas; // enter comment here
- Float_t fEtNonRemovedElectrons; // enter comment here
- Float_t fEtNonRemovedPositrons; // enter comment here
- Float_t fEtNonRemovedMuMinus; // enter comment here
- Float_t fEtNonRemovedMuPlus; // enter comment here
- Float_t fEtNonRemovedGammas; // enter comment here
- Float_t fEtNonRemovedGammasFromPi0; // enter comment here
- Float_t fEtNonRemovedNeutrons; // enter comment here
- Float_t fEtNonRemovedAntiNeutrons; // enter comment here
-
- Float_t fEtRemovedGammas; // enter comment here
- Float_t fEtRemovedNeutrons; // enter comment here
- Float_t fEtRemovedAntiNeutrons; // enter comment here
-
- Int_t fMultNonRemovedProtons; // enter comment here
- Int_t fMultNonRemovedAntiProtons; // enter comment here
- Int_t fMultNonRemovedPiPlus; // enter comment here
- Int_t fMultNonRemovedPiMinus; // enter comment here
- Int_t fMultNonRemovedKaonPlus; // enter comment here
- Int_t fMultNonRemovedKaonMinus; // enter comment here
- Int_t fMultNonRemovedK0s; // enter comment here
- Int_t fMultNonRemovedLambdas; // enter comment here
- Int_t fMultNonRemovedElectrons; // enter comment here
- Int_t fMultNonRemovedPositrons; // enter comment here
- Int_t fMultNonRemovedMuMinus; // enter comment here
- Int_t fMultNonRemovedMuPlus; // enter comment here
- Int_t fMultNonRemovedGammas; // enter comment here
- Int_t fMultNonRemovedNeutrons; // enter comment here
- Int_t fMultNonRemovedAntiNeutrons; // enter comment here
-
- Int_t fMultRemovedGammas; // enter comment here
- Int_t fMultRemovedNeutrons; // enter comment here
- Int_t fMultRemovedAntiNeutrons; // enter comment here
-
- Int_t fTrackMultInAcc; // enter comment here
-
-
- TH2F *fHistDxDzNonRemovedCharged; // enter comment here
- TH2F *fHistDxDzRemovedCharged; // enter comment here
- TH2F *fHistDxDzNonRemovedNeutral; // enter comment here
- TH2F *fHistDxDzRemovedNeutral; // enter comment here
-
- TH1F *fHistPiPlusMult; // enter comment here
- TH1F *fHistPiMinusMult; // enter comment here
- TH1F *fHistPiZeroMult; // enter comment here
-
- TH1F *fHistPiPlusMultAcc; // enter comment here
- TH1F *fHistPiMinusMultAcc; // enter comment here
- TH1F *fHistPiZeroMultAcc; // enter comment here
-
- Int_t fPiPlusMult; // enter comment here
- Int_t fPiMinusMult; // enter comment here
- Int_t fPiZeroMult; // enter comment here
-
- Int_t fPiPlusMultAcc; // enter comment here
- Int_t fPiMinusMultAcc; // enter comment here
- Int_t fPiZeroMultAcc; // enter comment here
-
-
- Int_t fNeutralRemoved; // number of neutral particles that where removed by track matching
- Int_t fChargedRemoved; // number of charged particles that where removed by track matching
- Int_t fChargedNotRemoved; // number of charged particles that were not removed
- Int_t fNeutralNotRemoved; // number of neutral particles that were not removed
-
- Double_t fEnergyNeutralRemoved; // energy of neutral particles that where removed by track matching
- Double_t fEnergyChargedRemoved; // energy of charged particles that where removed by track matching
- Double_t fEnergyChargedNotRemoved; // energy of charged particles that were not removed
- Double_t fEnergyNeutralNotRemoved; // energy of neutral particles that were not removed
-
-
- private:
+ TH3F *fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
+ TH3F *fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
+ TH3F *fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
+ TH3F *fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
+
+ TH2F *fHistRemovedOrNot; // If charged/neutral particles were removed or not
+
+ TH2F *fHistEtNonRemovedProtons; // enter comment here
+ TH2F *fHistEtNonRemovedAntiProtons; // enter comment here
+ TH2F *fHistEtNonRemovedPiPlus; // enter comment here
+ TH2F *fHistEtNonRemovedPiMinus; // enter comment here
+ TH2F *fHistEtNonRemovedKaonPlus; // enter comment here
+ TH2F *fHistEtNonRemovedKaonMinus; // enter comment here
+ TH2F *fHistEtNonRemovedK0s; // enter comment here
+ TH2F *fHistEtNonRemovedK0L; // enter comment here
+ TH2F *fHistEtNonRemovedLambdas; // enter comment here
+ TH2F *fHistEtNonRemovedElectrons; // enter comment here
+ TH2F *fHistEtNonRemovedPositrons; // enter comment here
+ TH2F *fHistEtNonRemovedMuPlus; // enter comment here
+ TH2F *fHistEtNonRemovedMuMinus; // enter comment here
+ TH2F *fHistEtNonRemovedNeutrons; // enter comment here
+ TH2F *fHistEtNonRemovedAntiNeutrons; // enter comment here
+ TH2F *fHistEtNonRemovedGammas; // enter comment here
+ TH2F *fHistEtNonRemovedGammasFromPi0; // enter comment here
+
+ TH2F *fHistEtRemovedGammas; // enter comment here
+ TH2F *fHistEtRemovedNeutrons; // enter comment here
+ TH2F *fHistEtRemovedAntiNeutrons; // enter comment here
+
+ TH2F *fHistEtRemovedCharged; // enter comment here
+ TH2F *fHistEtRemovedNeutrals; // enter comment here
+
+ TH2F *fHistEtNonRemovedCharged; // enter comment here
+ TH2F *fHistEtNonRemovedNeutrals; // enter comment here
+
+ TH2F *fHistMultNonRemovedProtons; // enter comment here
+ TH2F *fHistMultNonRemovedAntiProtons; // enter comment here
+ TH2F *fHistMultNonRemovedPiPlus; // enter comment here
+ TH2F *fHistMultNonRemovedPiMinus; // enter comment here
+ TH2F *fHistMultNonRemovedKaonPlus; // enter comment here
+ TH2F *fHistMultNonRemovedKaonMinus; // enter comment here
+ TH2F *fHistMultNonRemovedK0s; // enter comment here
+ TH2F *fHistMultNonRemovedK0L; // enter comment here
+ TH2F *fHistMultNonRemovedLambdas; // enter comment here
+ TH2F *fHistMultNonRemovedElectrons; // enter comment here
+ TH2F *fHistMultNonRemovedPositrons; // enter comment here
+ TH2F *fHistMultNonRemovedMuPlus; // enter comment here
+ TH2F *fHistMultNonRemovedMuMinus; // enter comment here
+ TH2F *fHistMultNonRemovedNeutrons; // enter comment here
+ TH2F *fHistMultNonRemovedAntiNeutrons; // enter comment here
+ TH2F *fHistMultNonRemovedGammas; // enter comment here
+
+ TH2F *fHistMultRemovedGammas; // enter comment here
+ TH2F *fHistMultRemovedNeutrons; // enter comment here
+ TH2F *fHistMultRemovedAntiNeutrons; // enter comment here
+
+ TH2F *fHistMultRemovedCharged; // enter comment here
+ TH2F *fHistMultRemovedNeutrals; // enter comment here
+
+ TH2F *fHistMultNonRemovedCharged; // enter comment here
+ TH2F *fHistMultNonRemovedNeutrals; // enter comment here
+
+ TH2F *fHistTrackMultvsNonRemovedCharged; // enter comment here
+ TH2F *fHistTrackMultvsNonRemovedNeutral; // enter comment here
+ TH2F *fHistTrackMultvsRemovedGamma; // enter comment here
+
+ TH2F *fHistClusterMultvsNonRemovedCharged; // enter comment here
+ TH2F *fHistClusterMultvsNonRemovedNeutral; // enter comment here
+ TH2F *fHistClusterMultvsRemovedGamma; // enter comment here
+
+ TH2F *fHistMultvsNonRemovedChargedE; // enter comment here
+ TH2F *fHistMultvsNonRemovedNeutralE; // enter comment here
+ TH2F *fHistMultvsRemovedGammaE; // enter comment here
+
+ Float_t fEtNonRemovedProtons; // enter comment here
+ Float_t fEtNonRemovedAntiProtons; // enter comment here
+ Float_t fEtNonRemovedPiPlus; // enter comment here
+ Float_t fEtNonRemovedPiMinus; // enter comment here
+ Float_t fEtNonRemovedKaonPlus; // enter comment here
+ Float_t fEtNonRemovedKaonMinus; // enter comment here
+ Float_t fEtNonRemovedK0S; // enter comment here
+ Float_t fEtNonRemovedK0L; // enter comment here
+ Float_t fEtNonRemovedLambdas; // enter comment here
+ Float_t fEtNonRemovedElectrons; // enter comment here
+ Float_t fEtNonRemovedPositrons; // enter comment here
+ Float_t fEtNonRemovedMuMinus; // enter comment here
+ Float_t fEtNonRemovedMuPlus; // enter comment here
+ Float_t fEtNonRemovedGammas; // enter comment here
+ Float_t fEtNonRemovedGammasFromPi0; // enter comment here
+ Float_t fEtNonRemovedNeutrons; // enter comment here
+ Float_t fEtNonRemovedAntiNeutrons; // enter comment here
+
+ Float_t fEtRemovedProtons; // enter comment here
+ Float_t fEtRemovedAntiProtons; // enter comment here
+ Float_t fEtRemovedPiPlus; // enter comment here
+ Float_t fEtRemovedPiMinus; // enter comment here
+ Float_t fEtRemovedKaonPlus; // enter comment here
+ Float_t fEtRemovedKaonMinus; // enter comment here
+ Float_t fEtRemovedK0s; // enter comment here
+ Float_t fEtRemovedK0L; // enter comment here
+ Float_t fEtRemovedLambdas; // enter comment here
+ Float_t fEtRemovedElectrons; // enter comment here
+ Float_t fEtRemovedPositrons; // enter comment here
+ Float_t fEtRemovedMuMinus; // enter comment here
+ Float_t fEtRemovedMuPlus; // enter comment here
+
+ Float_t fEtRemovedGammasFromPi0; // enter comment here
+ Float_t fEtRemovedGammas; // enter comment here
+ Float_t fEtRemovedNeutrons; // enter comment here
+ Float_t fEtRemovedAntiNeutrons; // enter comment here
+
+ Int_t fMultNonRemovedProtons; // enter comment here
+ Int_t fMultNonRemovedAntiProtons; // enter comment here
+ Int_t fMultNonRemovedPiPlus; // enter comment here
+ Int_t fMultNonRemovedPiMinus; // enter comment here
+ Int_t fMultNonRemovedKaonPlus; // enter comment here
+ Int_t fMultNonRemovedKaonMinus; // enter comment here
+ Int_t fMultNonRemovedK0s; // enter comment here
+ Int_t fMultNonRemovedK0L; // enter comment here
+ Int_t fMultNonRemovedLambdas; // enter comment here
+ Int_t fMultNonRemovedElectrons; // enter comment here
+ Int_t fMultNonRemovedPositrons; // enter comment here
+ Int_t fMultNonRemovedMuMinus; // enter comment here
+ Int_t fMultNonRemovedMuPlus; // enter comment here
+ Int_t fMultNonRemovedGammas; // enter comment here
+ Int_t fMultNonRemovedNeutrons; // enter comment here
+ Int_t fMultNonRemovedAntiNeutrons; // enter comment here
+
+ Int_t fMultRemovedProtons; // enter comment here
+ Int_t fMultRemovedAntiProtons; // enter comment here
+ Int_t fMultRemovedPiPlus; // enter comment here
+ Int_t fMultRemovedPiMinus; // enter comment here
+ Int_t fMultRemovedKaonPlus; // enter comment here
+ Int_t fMultRemovedKaonMinus; // enter comment here
+ Int_t fMultRemovedK0s; // enter comment here
+ Int_t fMultRemovedK0L; // enter comment here
+
+ Int_t fMultRemovedLambdas; // enter comment here
+ Int_t fMultRemovedElectrons; // enter comment here
+ Int_t fMultRemovedPositrons; // enter comment here
+ Int_t fMultRemovedMuMinus; // enter comment here
+ Int_t fMultRemovedMuPlus; // enter comment here
+
+ Int_t fMultRemovedGammas; // enter comment here
+ Int_t fMultRemovedNeutrons; // enter comment here
+ Int_t fMultRemovedAntiNeutrons; // enter comment here
+
+ Int_t fTrackMultInAcc; // enter comment here
+
+
+ TH2F *fHistDxDzNonRemovedCharged; // enter comment here
+ TH2F *fHistDxDzRemovedCharged; // enter comment here
+ TH2F *fHistDxDzNonRemovedNeutral; // enter comment here
+ TH2F *fHistDxDzRemovedNeutral; // enter comment here
+
+ TH1F *fHistPiPlusMult; // enter comment here
+ TH1F *fHistPiMinusMult; // enter comment here
+ TH1F *fHistPiZeroMult; // enter comment here
+
+ TH1F *fHistPiPlusMultAcc; // enter comment here
+ TH1F *fHistPiMinusMultAcc; // enter comment here
+ TH1F *fHistPiZeroMultAcc; // enter comment here
+
+ Int_t fPiPlusMult; // enter comment here
+ Int_t fPiMinusMult; // enter comment here
+ Int_t fPiZeroMult; // enter comment here
+
+ Int_t fPiPlusMultAcc; // enter comment here
+ Int_t fPiMinusMultAcc; // enter comment here
+ Int_t fPiZeroMultAcc; // enter comment here
+
+
+ Int_t fNeutralRemoved; // number of neutral particles that where removed by track matching
+ Int_t fChargedRemoved; // number of charged particles that where removed by track matching
+ Int_t fChargedNotRemoved; // number of charged particles that were not removed
+ Int_t fNeutralNotRemoved; // number of neutral particles that were not removed
+
+ Double_t fEnergyNeutralRemoved; // energy of neutral particles that where removed by track matching
+ Double_t fEnergyChargedRemoved; // energy of charged particles that where removed by track matching
+ Double_t fEnergyChargedNotRemoved; // energy of charged particles that were not removed
+ Double_t fEnergyNeutralNotRemoved; // energy of neutral particles that were not removed
+
+ Int_t fNClusters; // Number of clusters in event
+
+ AliPHOSGeometry *fGeoUtils;
+
+ TH2I *fBadMapM2; // Bad map
+ TH2I *fBadMapM3; // Bad map
+ TH2I *fBadMapM4; // Bad map
+
+private:
//Declare it private to avoid compilation warning
AliAnalysisEtMonteCarlo & operator = (const AliAnalysisEtMonteCarlo & g) ;//cpy assignment
AliAnalysisEtMonteCarloEmcal::AliAnalysisEtMonteCarloEmcal()
{
fHistogramNameSuffix = TString("EmcalMC");
+
}
AliAnalysisEtMonteCarloEmcal::~AliAnalysisEtMonteCarloEmcal()
void AliAnalysisEtMonteCarloEmcal::Init()
{ // Init
AliAnalysisEtMonteCarlo::Init();
-
+ fClusterType = fCuts->GetEmcalClusterType();
fDetectorRadius = fCuts->GetGeometryEmcalDetectorRadius();
fEtaCutAcc = fCuts->GetGeometryEmcalEtaAccCut();
fPhiCutAccMax = fCuts->GetGeometryEmcalPhiAccMaxCut() * TMath::Pi()/180.;
//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
//_________________________________________________________________________
#include "AliAnalysisEtMonteCarloPhos.h"
+#include "AliAnalysisEtSelectorPhos.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDtrack.h"
+#include <iostream>
+#include "AliPHOSGeoUtils.h"
+#include "TFile.h"
+#include "TH2I.h"
+#include <AliPHOSGeometry.h>
using namespace std;
ClassImp(AliAnalysisEtMonteCarloPhos);
-AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos()
+AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos():AliAnalysisEtMonteCarlo()
+,fBadMapM2(0)
+,fBadMapM3(0)
+,fBadMapM4(0)
{
fHistogramNameSuffix = TString("PhosMC");
}
void AliAnalysisEtMonteCarloPhos::Init()
{ // Init
AliAnalysisEtMonteCarlo::Init();
+ fSelector = new AliAnalysisEtSelectorPhos(fCuts);
+ fSelector->Init(137366);
+ fClusterType = fCuts->GetPhosClusterType();
fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
fEtaCutAcc = fCuts->GetGeometryPhosEtaAccCut();
fPhiCutAccMax = fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180.;
fTrackDzCut = fCuts->GetPhosTrackDzCut();
fDetector = fCuts->GetDetectorPhos();
-
+
+ // ifstream f("badchannels.txt", ios::in);
+ TFile *f = TFile::Open("badchannels.root", "READ");
+
+ fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
+ fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
+ fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
+//
+// fGeoUtils = new AliPHOSGeoUtils("PHOS", "noCPV");
+ fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
+
}
+
+
+Bool_t AliAnalysisEtMonteCarloPhos::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
+{
+
+ Float_t gPos[3];
+ cluster.GetPosition(gPos);
+ Int_t relId[4];
+ TVector3 glVec(gPos);
+ fGeoUtils->GlobalPos2RelId(glVec, relId);
+
+ TVector3 locVec;
+ fGeoUtils->Global2Local(locVec, glVec, relId[0]);
+// std::cout << fGeoUtils << std::endl;
+ //std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
+ //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
+ for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
+ {
+ for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
+ {
+ if (relId[0] == 3)
+ {
+ if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 3;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
+
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 2)
+ {
+ if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 2;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
+
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 1)
+ {
+ if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 1;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
+
+ return kTRUE;
+ }
+ }
+ }
+
+ }
+ }
+
+ return kFALSE;
+}
+
+
#ifndef ALIANALYSISETMONTECARLOPHOS_H
#define ALIANALYSISETMONTECARLOPHOS_H
-//_________________________________________________________________________
+
+class AliPHOSGeoUtils;//_________________________________________________________________________
// Utility Class for transverse energy studies
// Base class for MC analysis, for PHOS
// - MC output
//_________________________________________________________________________
#include "AliAnalysisEtMonteCarlo.h"
-
+class TH2I;
class AliAnalysisEtMonteCarloPhos : public AliAnalysisEtMonteCarlo
{
virtual ~AliAnalysisEtMonteCarloPhos();
virtual void Init();
+protected:
+ virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
private:
-
+
+ TH2I *fBadMapM2; // Bad map
+ TH2I *fBadMapM3; // Bad map
+ TH2I *fBadMapM4; // Bad map
+
+ // Prohibited
+ AliAnalysisEtMonteCarloPhos & operator = (const AliAnalysisEtMonteCarloPhos&) ;//cpy assignment
+ AliAnalysisEtMonteCarloPhos(const AliAnalysisEtMonteCarloPhos&) ; // cpy ctor
ClassDef(AliAnalysisEtMonteCarloPhos, 1);
};
#include "AliESDpid.h"
#include <iostream>
#include "TH2F.h"
+#include "TH2I.h"
+#include "TH1I.h"
+#include "TFile.h"
#include "AliAnalysisHadEtCorrections.h"
+#include "AliAnalysisEtSelector.h"
#include "AliLog.h"
#include "AliCentrality.h"
-
-
+#include "AliPHOSGeoUtils.h"
+#include "AliPHOSGeometry.h"
using namespace std;
AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
- AliAnalysisEt()
- ,fCorrections(0)
- ,fPidCut(0)
- ,fClusterType(0)
- ,fHistChargedPionEnergyDeposit(0)
- ,fHistProtonEnergyDeposit(0)
- ,fHistAntiProtonEnergyDeposit(0)
- ,fHistChargedKaonEnergyDeposit(0)
- ,fHistMuonEnergyDeposit(0)
- ,fHistRemovedEnergy(0)
- ,fGeomCorrection(1.0)
- ,fEMinCorrection(1.0)
+ AliAnalysisEt()
+ ,fCorrections(0)
+ ,fPidCut(0)
+ ,fHistChargedPionEnergyDeposit(0)
+ ,fHistProtonEnergyDeposit(0)
+ ,fHistAntiProtonEnergyDeposit(0)
+ ,fHistChargedKaonEnergyDeposit(0)
+ ,fHistMuonEnergyDeposit(0)
+ ,fHistRemovedEnergy(0)
+ ,fGeomCorrection(1.0)
+ ,fEMinCorrection(1.0/0.89)
+ ,fRecEffCorrection(1.0)
+ ,fGeoUtils(0)
+ ,fBadMapM2(0)
+ ,fBadMapM3(0)
+ ,fBadMapM4(0)
+ ,fClusterPosition(0)
+ ,fHistChargedEnergyRemoved(0)
+ ,fHistNeutralEnergyRemoved(0)
+ ,fHistGammaEnergyAdded(0)
{
}
AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
-{//destructor
+{ //destructor
delete fCorrections;
- delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
- delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
- delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
- delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
+ delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
+ delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
+ delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
+ delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
delete fHistRemovedEnergy; // removed energy
AliFatal("ERROR: ESD Event does not exist");
return 0;
}
+ fSelector->SetEvent(event);
+ if (!fMatrixInitialized)
+ {
+ for (Int_t mod=0; mod<5; mod++) {
+ if (!event->GetPHOSMatrix(mod)) continue;
+ fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
+// std::cout << event->GetPHOSMatrix(mod) << std::endl;
+ Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
+ }
+ fMatrixInitialized = kTRUE;
+ }
Int_t cent = -1;
if (fCentrality)
{
cent = fCentrality->GetCentralityClass10("V0M");
- fCentClass = fCentrality->GetCentralityClass10("V0M");
+ fCentClass = fCentrality->GetCentralityClass10("V0M");
}
- Double_t protonMass = fgProtonMass;
-
- //for PID
- AliESDpid pID;
- pID.MakePID(event);
- TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(event);
-
- Int_t nGoodTracks = list->GetEntries();
- // printf("nGoodTracks %d nCaloClusters %d\n", nGoodTracks, event->GetNumberOfCaloClusters());
-
- for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
- {
- AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
- if (!track)
- {
- AliError(Form("ERROR: Could not get track %d", iTrack));
- continue;
- }
-
- fMultiplicity++;
-
-
- Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
- nSigmaPion = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kPion));
- nSigmaProton = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kProton));
- nSigmaKaon = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kKaon));
- nSigmaElectron = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kElectron));
- /*
- bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
- bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
- bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
- bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
- */
-
- Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
- Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
-
- Float_t massPart = 0;
-
- const Double_t *pidWeights = track->PID();
- Int_t maxpid = -1;
- Double_t maxpidweight = 0;
-
- if (pidWeights)
- {
- for (Int_t p =0; p < AliPID::kSPECIES; p++)
- {
- if (pidWeights[p] > maxpidweight)
- {
- maxpidweight = pidWeights[p];
- maxpid = p;
- }
- }
- if (maxpid == AliPID::kProton)
- {
- //by definition of ET
- massPart = -protonMass*track->Charge();
- }
-
- }
-
- Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
- if(fMakeSparse){
- fSparseTracks[0] = maxpid;
- fSparseTracks[1] = track->Charge();
- fSparseTracks[2] = track->M();
- fSparseTracks[3] = et;
- fSparseTracks[4] = track->Pt();
- fSparseTracks[5] = track->Eta();
- fSparseTracks[6] = cent;
- fSparseHistTracks->Fill(fSparseTracks);
- }
- //printf("Rec track: iTrack %03d eta %4.3f phi %4.3f nITSCl %d nTPCCl %d\n", iTrack, track->Eta(), track->Phi(), nItsClusters, nTPCClusters); // tmp/debug printout
-
- if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
- {
- fTotChargedEt += et;
- fChargedMultiplicity++;
- if (maxpid != -1)
- {
- if (maxpid == AliPID::kProton)
- {
- fProtonEt += et;
- }
- if (maxpid == AliPID::kPion)
- {
- fPionEt += et;
- }
- if (maxpid == AliPID::kKaon)
- {
- fChargedKaonEt += et;
- }
- if (maxpid == AliPID::kMuon)
- {
- fMuonEt += et;
- }
- if (maxpid == AliPID::kElectron)
- {
- fElectronEt += et;
- }
- }
-
- if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
- {
- fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
- if (maxpid != -1)
- {
- if (maxpid == AliPID::kProton)
- {
- fProtonEtAcc += et;
- }
- if (maxpid == AliPID::kPion)
- {
- fPionEtAcc += et;
- }
- if (maxpid == AliPID::kKaon)
- {
- fChargedKaonEtAcc += et;
- }
- if (maxpid == AliPID::kMuon)
- {
- fMuonEtAcc += et;
- }
- if (maxpid == AliPID::kElectron)
- {
- fElectronEtAcc += et;
- }
- }
-
- }
- }
-
- if (TrackHitsCalorimeter(track, event->GetMagneticField()))
- {
- Double_t phi = track->Phi();
- Double_t pt = track->Pt();
- // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
- if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
- else fHistPhivsPtNeg->Fill(phi, pt);
- }
- }
+ //Double_t protonMass = fgProtonMass;
for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
{
AliError(Form("ERROR: Could not get cluster %d", iCluster));
continue;
}
- if (cluster->GetType() != fClusterType) continue;
-
- //if(cluster->GetTracksMatched() > 0)
-// printf("Rec Cluster: iCluster %03d E %4.3f type %.qd NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout
-
-
- if (cluster->E() < fClusterEnergyCut) continue;
+ int x = 0;
+ fCutFlow->Fill(x++);
+ if(cluster->IsEMCAL()) continue;
+ fCutFlow->Fill(x++);
+ if(!fSelector->CutMinEnergy(*cluster)) continue;
+ fCutFlow->Fill(x++);
+ if (!fSelector->CutDistanceToBadChannel(*cluster)) continue;
+ fCutFlow->Fill(x++);
Float_t pos[3];
distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
}
- if(fMakeSparse){
- fSparseClusters[0] = 0;
- fSparseClusters[1] = 0;
- fSparseClusters[2] = 0;
- fSparseClusters[6] = 0;
- fSparseClusters[7] = 0;
- fSparseClusters[8] = 0;
- fSparseClusters[9] = cent;
- fSparseClusters[10] = 0;
- }
-
- if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex > -1)
- {
- AliVTrack *tmptrack = event->GetTrack(trackMatchedIndex);
- if (!tmptrack)
- {
- AliError("Error: track does not exist");
- return -1;
- }
- const Double_t *pidWeights = tmptrack->PID();
-
- Double_t maxpidweight = 0;
- Int_t maxpid = 0;
- Double_t massPart = 0;
- if (pidWeights)
- {
- for (Int_t p =0; p < AliPID::kSPECIES; p++)
- {
- if (pidWeights[p] > maxpidweight)
- {
- maxpidweight = pidWeights[p];
- maxpid = p;
- }
- }
- if (maxpid == AliPID::kProton)
- {
- //by definition of ET
- massPart = -protonMass*tmptrack->Charge();
- }
- }
- if(fMakeSparse){
- fSparseClusters[0] = maxpid;
- fSparseClusters[1] = tmptrack->Charge();
- fSparseClusters[2] = tmptrack->M();
- fSparseClusters[6] = tmptrack->E() * TMath::Sin(tmptrack->Theta()) + massPart;;
- fSparseClusters[7] = tmptrack->Pt();
- fSparseClusters[8] = tmptrack->Eta();
- }
- }
-
-
- if(fMakeSparse){fSparseClusters[10] = distance;}
-
- fHistTMDeltaR->Fill(distance);
- fHistTMDxDz->Fill(cluster->GetTrackDx(), cluster->GetTrackDz());
-
-// Float_t clusteret = cluster->E() * TMath::Sin(cp.Theta());
-
Bool_t matched = false;
- if (cluster->IsEMCAL()) matched = distance < fTrackDistanceCut;
- else matched = (TMath::Abs(cluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(cluster->GetTrackDz()) < fTrackDzCut);
+ Double_t r = 0;
+
+ matched = !fSelector->CutTrackMatching(*cluster, r);
if (matched)
{
+
if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
{
AliVTrack *track = event->GetTrack(trackMatchedIndex);
} // distance
else
{
- if(fMakeSparse){
- fSparseClusters[0] = AliPID::kPhoton;
- fSparseClusters[1] = 0;
- }
+ fCutFlow->Fill(x++);
+ //std::cout << x++ << std::endl;
+ fSparseClusters[0] = AliPID::kPhoton;
+ fSparseClusters[1] = 0;
- if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
- if (cluster->E() < fClusterEnergyCut) continue;
+ //if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
+ //if (cluster->E() < fClusterEnergyCut) continue;
cluster->GetPosition(pos);
- // TODO: replace with TVector3, too lazy now...
+ TVector3 p2(pos);
- float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
+ fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
- float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
- // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
- fTotNeutralEt += cluster->E() * TMath::Sin(theta);
+ fTotNeutralEt += CalculateTransverseEnergy(cluster);
fNeutralMultiplicity++;
}
+ fMultiplicity++;
+ }
- cluster->GetPosition(pos);
-
- // TODO: replace with TVector3
+ fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
+ fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
+ fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
+ fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
- float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
- float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
- float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
- if(fMakeSparse){
- fSparseClusters[3] = cluster->E() * TMath::Sin(theta);
- fSparseClusters[4] = cluster->E();
- fSparseClusters[5] = eta;
-
- fSparseHistClusters->Fill(fSparseClusters);
- }
+ fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
+ fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
- fMultiplicity++;
- }
Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity);
fHistRemovedEnergy->Fill(removedEnergy);
-// std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
- fTotNeutralEt = fGeomCorrection * fEMinCorrection * fTotNeutralEt - removedEnergy;
- fTotNeutralEtAcc = fTotNeutralEt/fGeomCorrection;
+
+ fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
+ fTotNeutralEtAcc = fTotNeutralEt;
fTotEt = fTotChargedEt + fTotNeutralEt;
fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
- if(fMakeSparse){
- fSparseEt[0] = fTotEt;
- fSparseEt[1] = fTotNeutralEt;
- fSparseEt[2] = fTotChargedEtAcc;
- fSparseEt[3] = fMultiplicity;
- fSparseEt[4] = fNeutralMultiplicity;
- fSparseEt[5] = fChargedMultiplicity;
- fSparseEt[6] = cent;
- }
+ if(fMakeSparse) {
+ fSparseEt[0] = fTotEt;
+ fSparseEt[1] = fTotNeutralEt;
+ fSparseEt[2] = fTotChargedEtAcc;
+ fSparseEt[3] = fMultiplicity;
+ fSparseEt[4] = fNeutralMultiplicity;
+ fSparseEt[5] = fChargedMultiplicity;
+ fSparseEt[6] = cent;
+ }
// Fill the histograms...
FillHistograms();
}
bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
-{ // check vertex
+{ // check vertex
Float_t bxy = 999.;
Float_t bz = 999.;
}
void AliAnalysisEtReconstructed::Init()
-{ // Init
+{ // Init
AliAnalysisEt::Init();
fPidCut = fCuts->GetReconstructedPidCut();
TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
if (!fCorrections) {
cout<<"Warning! You have not set corrections. Your code will crash. You have to set the corrections."<<endl;
}
+ //fGeoUtils = new AliPHOSGeoUtils("PHOS", "noCPV");
+ fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
+ // ifstream f("badchannels.txt", ios::in);
+ TFile *f = TFile::Open("badchannels.root", "READ");
+
+ fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
+ fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
+ fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
}
bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
-{ // propagate track to detector radius
+{ // propagate track to detector radius
if (!track) {
cout<<"Warning: track empty"<<endl;
}
void AliAnalysisEtReconstructed::FillOutputList(TList* list)
-{ // add some extra histograms to the ones from base class
+{ // add some extra histograms to the ones from base class
AliAnalysisEt::FillOutputList(list);
list->Add(fHistChargedPionEnergyDeposit);
list->Add(fHistMuonEnergyDeposit);
list->Add(fHistRemovedEnergy);
+ list->Add(fClusterPosition);
+
+ list->Add(fHistChargedEnergyRemoved);
+ list->Add(fHistNeutralEnergyRemoved);
+ list->Add(fHistGammaEnergyAdded);
}
void AliAnalysisEtReconstructed::CreateHistograms()
-{ // add some extra histograms to the ones from base class
+{ // add some extra histograms to the ones from base class
AliAnalysisEt::CreateHistograms();
Int_t nbinsEt = 1000;
//fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
//fHistMuonEnergyDeposit->SetYTitle("Energy of track");
+ histname = "fClusterPosition" + fHistogramNameSuffix;
+ fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",1000, -2.0, -.5, 1000, -.13 , 0.13);
+ fClusterPosition->SetXTitle("Energy deposited in calorimeter");
+ fClusterPosition->SetYTitle("Energy of track");
+
+
+ histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
+ fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
+
+ histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
+ fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
+
+ histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
+ fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
}
AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
Int_t *trkMatchId,
const AliESDEvent *event)
-{ // calculate distance between cluster and closest track
+{ // calculate distance between cluster and closest track
Double_t trkPos[3] = {0,0,0};
return distance;
}
+
+Bool_t AliAnalysisEtReconstructed::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
+{
+
+ Float_t gPos[3];
+ cluster.GetPosition(gPos);
+ Int_t relId[4];
+ TVector3 glVec(gPos);
+ fGeoUtils->GlobalPos2RelId(glVec, relId);
+
+ TVector3 locVec;
+ fGeoUtils->Global2Local(locVec, glVec, relId[0]);
+// std::cout << fGeoUtils << std::endl;
+ //std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
+ //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
+ for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
+ {
+ for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
+ {
+ if (relId[0] == 3)
+ {
+ if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 3;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 2)
+ {
+ if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 2;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 1)
+ {
+ if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 1;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kTRUE;
+ }
+ }
+ }
+
+ }
+ }
+
+ return kFALSE;
+
+
+}
+
+
+
+
+/*
+Bool_t AliAnalysisEtReconstructed::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
+{
+
+ Float_t gPos[3];
+
+ cluster.GetPosition(gPos);
+ Int_t relId[4];
+ TVector3 glVec(gPos);
+ fGeoUtils->GlobalPos2RelId(glVec, relId);
+ TVector3 locVec;
+ fGeoUtils->Global2Local(locVec, glVec, relId[0]);
+
+ std::vector<Int_t>::const_iterator badIt;
+
+ for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
+ {
+ for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
+ {
+ if (relId[0] == 3)
+ {
+ if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
+ {
+
+ Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 2)
+ {
+ if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
+ {
+
+ Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+ return kTRUE;
+ }
+ }
+ }
+ if (relId[0] == 1)
+ {
+ if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
+ {
+
+ Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+ return kTRUE;
+ }
+ }
+ }
+ }
+ }
+
+ return kFALSE;
+}
+*/
#ifndef ALIANALYSISETRECONSTRUCTED_H
#define ALIANALYSISETRECONSTRUCTED_H
-//_________________________________________________________________________
+
+class TH2F;
+class TH2D;
+class AliPHOSGeometry;
+class AliPHOSGeoUtils;//_________________________________________________________________________
// Utility Class for transverse energy studies
// Base class for ESD analysis
// - reconstruction output
#include "AliAnalysisEt.h"
+class TH2I;
class AliVParticle;
class AliESDEvent;
class AliAnalysisHadEtCorrections;
{
public:
-
+
AliAnalysisEtReconstructed();
virtual ~AliAnalysisEtReconstructed();
/** Fill the objects you want to output, classes which add new histograms should overload this. */
virtual void FillOutputList(TList *list);
- void SetCorrections(AliAnalysisHadEtCorrections *corr){fCorrections = corr;}
+ void SetCorrections(AliAnalysisHadEtCorrections *corr) {
+ fCorrections = corr;
+ }
/** Create the histograms, must be overloaded if you want to add your own */
virtual void CreateHistograms();
+ void SetEMinCorrection(const Double_t factor) { fEMinCorrection = factor; }
+
protected:
bool CheckGoodVertex(AliVParticle *track);
virtual bool TrackHitsCalorimeter(AliVParticle *track, Double_t magField);
+ virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
+
+ //Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
AliAnalysisHadEtCorrections *fCorrections;//corrections needed for hadronic et
Double_t fPidCut; // cut on the pid probability
-
- Char_t fClusterType; // selection on cluster type
-
- TH2F *fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
- TH2F *fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
- TH2F *fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
- TH2F *fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
+
+ TH2F *fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
+ TH2F *fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
+ TH2F *fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
+ TH2F *fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
TH2F *fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
TH1F *fHistRemovedEnergy; // removed energy
-
+
Double_t fGeomCorrection; // geometry correction
Double_t fEMinCorrection; // Emin correction
- private:
+ Double_t fRecEffCorrection; // Eff correction
+ Double_t CalcTrackClusterDistance(const Float_t pos[3],Int_t *trkMatchId, const AliESDEvent *event);
+
+// AliPHOSGeoUtils *fGeoUtils;
+ AliPHOSGeometry *fGeoUtils;
+ TH2I *fBadMapM2; // Bad map
+ TH2I *fBadMapM3; // Bad map
+ TH2I *fBadMapM4; // Bad map
+
+
+ TH2D *fClusterPosition; // Position of clusters
+ TH2D *fHistChargedEnergyRemoved;
+ TH2D *fHistNeutralEnergyRemoved;
+ TH2D *fHistGammaEnergyAdded;
+
+
+private:
+
AliAnalysisEtReconstructed(const AliAnalysisEtReconstructed& g);
AliAnalysisEtReconstructed & operator=(const AliAnalysisEtReconstructed&);
-
- Double_t CalcTrackClusterDistance(const Float_t pos[3],Int_t *trkMatchId, const AliESDEvent *event);
+
+
ClassDef(AliAnalysisEtReconstructed, 1);
};
fClusterEnergyCut = fCuts->GetReconstructedEmcalClusterEnergyCut();
fSingleCellEnergyCut = fCuts->GetReconstructedEmcalSingleCellEnergyCut();
- fClusterType = fCuts->GetReconstructedEmcalClusterType();
+ fClusterType = fCuts->GetEmcalClusterType();
fTrackDistanceCut = fCuts->GetEmcalTrackDistanceCut();
fDetector = fCuts->GetDetectorEmcal();
//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
//_________________________________________________________________________
#include "AliAnalysisEtReconstructedPhos.h"
+#include "AliAnalysisEtSelectorPhos.h"
#include "AliAnalysisEtCuts.h"
#include "AliESDtrack.h"
+#include <fstream>
+#include <iostream>
+#include <AliPHOSGeoUtils.h>
+#include "AliPHOSGeometry.h"
+#include <vector>
+#include "TH2I.h"
+
using namespace std;
ClassImp(AliAnalysisEtReconstructedPhos);
const Double_t kMEANNEUTRAL = 0.434;
const Double_t kMEANGAMMA = 0.374;
-// Best case (pions and K0s):
+// Best case (pions and K0s):
const Double_t kMEANCHARGED = 0.304;
const Double_t kMEANNEUTRAL = 0.3356;
const Double_t kMEANGAMMA = 0.374;
*/
// Simulated case:
-const Double_t kMEANCHARGED = 0.307;
-const Double_t kMEANNEUTRAL = 0.407;
-const Double_t kMEANGAMMA = 0.374;
+//corrEnergy =cluster->E()/(0.51 + 0.02*cluster->E());
+const Double_t kMEANCHARGED = 0.307/(0.51 + 0.02*0.307);
+const Double_t kMEANNEUTRAL = 0.407/(0.51 + 0.02*0.407);
+const Double_t kMEANGAMMA = 0.374/(0.51 + 0.02*0.374);
AliAnalysisEtReconstructedPhos::AliAnalysisEtReconstructedPhos() :
-AliAnalysisEtReconstructed()
+ AliAnalysisEtReconstructed()
{
- fHistogramNameSuffix = TString("PhosRec");
+ fHistogramNameSuffix = TString("PhosRec");
+
+ fChargedContributionCorrectionParameters[0] = -0.017;
+ fChargedContributionCorrectionParameters[1] = 0.065;
+
+ fNeutralContributionCorrectionParameters[0] = -0.002;
+ fNeutralContributionCorrectionParameters[1] = 0.017;
+ fNeutralContributionCorrectionParameters[2] = -3.6e-5;
+
+ fRemovedGammaContributionCorrectionParameters[0] = 0.001;
+ fRemovedGammaContributionCorrectionParameters[1] = 0.37e-5;
+ fRemovedGammaContributionCorrectionParameters[2] = 0.0003;
+
}
-AliAnalysisEtReconstructedPhos::~AliAnalysisEtReconstructedPhos()
+AliAnalysisEtReconstructedPhos::~AliAnalysisEtReconstructedPhos()
{
}
void AliAnalysisEtReconstructedPhos::Init()
{ // Init
- AliAnalysisEtReconstructed::Init();
+ AliAnalysisEtReconstructed::Init();
+
+ fSelector = new AliAnalysisEtSelectorPhos(fCuts);
- fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
- fEtaCutAcc = fCuts->GetGeometryPhosEtaAccCut();
- fPhiCutAccMax = fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180.;
- fPhiCutAccMin = fCuts->GetGeometryPhosPhiAccMinCut() * TMath::Pi()/180.;
- fClusterEnergyCut = fCuts->GetReconstructedPhosClusterEnergyCut();
- fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();
-
- fClusterType = fCuts->GetReconstructedPhosClusterType();
- fTrackDistanceCut = fCuts->GetPhosTrackDistanceCut();
- fTrackDxCut = fCuts->GetPhosTrackDxCut();
- fTrackDzCut = fCuts->GetPhosTrackDzCut();
-
- fDetector = fCuts->GetDetectorPhos();
-
- fGeomCorrection = 1.0/0.036;
+ fSelector->Init(137366);
+
+ fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
+ fEtaCutAcc = fCuts->GetGeometryPhosEtaAccCut();
+ fPhiCutAccMax = fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180.;
+ fPhiCutAccMin = fCuts->GetGeometryPhosPhiAccMinCut() * TMath::Pi()/180.;
+ fClusterEnergyCut = fCuts->GetReconstructedPhosClusterEnergyCut();
+ fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();
+
+ fClusterType = fCuts->GetPhosClusterType();
+ fTrackDistanceCut = fCuts->GetPhosTrackDistanceCut();
+ fTrackDxCut = fCuts->GetPhosTrackDxCut();
+ fTrackDzCut = fCuts->GetPhosTrackDzCut();
+
+ fDetector = fCuts->GetDetectorPhos();
- fEMinCorrection = 1.0;
}
bool AliAnalysisEtReconstructedPhos::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
{
- return AliAnalysisEtReconstructed::TrackHitsCalorimeter(track, magField);
+ return AliAnalysisEtReconstructed::TrackHitsCalorimeter(track, magField);
}
Double_t AliAnalysisEtReconstructedPhos::GetChargedContribution(Int_t clusterMult)
{ // Charged contrib
- if(clusterMult > 0)
- {
- Double_t nPart = 0.067 + 0.137*clusterMult;
-
- Double_t contr = nPart*kMEANCHARGED;
-
- return contr;
- }
- return 0;
-
+ if (clusterMult > 0)
+ {
+ Double_t nPart = fChargedContributionCorrectionParameters[0] + fChargedContributionCorrectionParameters[1]*clusterMult;
+
+ Double_t contr = nPart*kMEANCHARGED;
+
+ return contr;
+ }
+ return 0;
+
}
Double_t AliAnalysisEtReconstructedPhos::GetNeutralContribution(Int_t clusterMult)
{ // Neutral contrib
- if(clusterMult > 0)
- {
- Double_t nPart = 0.012 + 0.024*clusterMult - 0.00006*clusterMult*clusterMult;
-
- Double_t contr = nPart*kMEANNEUTRAL;
-
- return contr;
- }
- return 0;
+ if (clusterMult > 0)
+ {
+ Double_t nPart = fNeutralContributionCorrectionParameters[0] + fNeutralContributionCorrectionParameters[1]*clusterMult + fNeutralContributionCorrectionParameters[2]*clusterMult*clusterMult;
+
+ Double_t contr = nPart*kMEANNEUTRAL;
+
+ return contr;
+ }
+ return 0;
}
Double_t AliAnalysisEtReconstructedPhos::GetGammaContribution(Int_t clusterMult)
{ // Gamma contrib
- if(clusterMult > 0)
- {
- Double_t nPart = -0.008 + 0.0057*clusterMult + 0.0002*clusterMult*clusterMult;
-
- Double_t contr = nPart*kMEANGAMMA;
-
- return contr;
- }
- return 0;
-}
+ if (clusterMult > 0)
+ {
+ Double_t nPart = fRemovedGammaContributionCorrectionParameters[0] + fRemovedGammaContributionCorrectionParameters[1]*clusterMult + fRemovedGammaContributionCorrectionParameters[2]*clusterMult*clusterMult;
+ Double_t contr = nPart*kMEANGAMMA;
+ return contr;
+ }
+ return 0;
+}
{
public:
-
+
AliAnalysisEtReconstructedPhos();
virtual ~AliAnalysisEtReconstructedPhos();
virtual void Init();
-
+
virtual Double_t GetChargedContribution(Int_t clusterMultiplicity);
-
+
virtual Double_t GetNeutralContribution(Int_t clusterMultiplicity);
-
+
virtual Double_t GetGammaContribution(Int_t clusterMultiplicity);
+
+ void SetChargedContributionParameters(Double_t par[2])
+ {
+ fChargedContributionCorrectionParameters[0] = par[0];
+ fChargedContributionCorrectionParameters[1] = par[1];
+ }
+ void SetNeutralContributionParameters(Double_t par[3])
+ {
+ fNeutralContributionCorrectionParameters[0] = par[0];
+ fNeutralContributionCorrectionParameters[1] = par[1];
+ fNeutralContributionCorrectionParameters[2] = par[2];
+ }
+ void SetRemovedGammaContributionParameters(Double_t par[3])
+ {
+ fRemovedGammaContributionCorrectionParameters[0] = par[0];
+ fRemovedGammaContributionCorrectionParameters[1] = par[1];
+ fRemovedGammaContributionCorrectionParameters[2] = par[2];
+ }
+
+protected:
+
+ virtual bool TrackHitsCalorimeter(AliVParticle *track, Double_t magField);
- protected:
-
- virtual bool TrackHitsCalorimeter(AliVParticle *track, Double_t magField);
- private:
+private:
+
+ Double_t fChargedContributionCorrectionParameters[2]; // Parametrization of the charged contribution as function of cluster multiplicity
+ Double_t fNeutralContributionCorrectionParameters[3]; // Parametrization of the neutral contribution as function of cluster multiplicity
+ Double_t fRemovedGammaContributionCorrectionParameters[3]; // Parametrization of the negative contribution from removed gammas as function of cluster multiplicity
+
- ClassDef(AliAnalysisEtReconstructedPhos, 1);
+ ClassDef(AliAnalysisEtReconstructedPhos, 1);
};
#endif // ALIANALYSISETRECONSTRUCTEDPHOS_H
--- /dev/null
+#include "AliAnalysisEtSelector.h"
+#include "AliAnalysisEtCuts.h"
+#include <AliESDCaloCluster.h>
+
+ClassImp(AliAnalysisEtSelector)
+
+AliAnalysisEtSelector::AliAnalysisEtSelector(AliAnalysisEtCuts *cuts) :
+fEvent(0)
+,fCuts(cuts)
+,fRunNumber(0)
+{
+
+}
+
+AliAnalysisEtSelector::~AliAnalysisEtSelector()
+{
+
+}
+
--- /dev/null
+#ifndef ALIANALYSISETSELECTOR_H
+#define ALIANALYSISETSELECTOR_H
+#include <Rtypes.h>
+
+class AliVEvent;
+class AliESDCaloCluster;
+class TRefArray;
+class AliAnalysisEtCuts;
+class AliAnalysisEtSelector
+{
+
+public:
+
+ // Constructor takes cuts object
+ AliAnalysisEtSelector(AliAnalysisEtCuts *cuts);
+
+ // Destructor
+ virtual ~AliAnalysisEtSelector();
+
+ // Set the current event
+ void SetEvent(const AliVEvent *event) { fEvent = event; }
+
+ // Init
+ virtual Int_t Init(Int_t runNumber) { fRunNumber = runNumber; return 0; }
+
+ // Return CaloClusters for the detector
+ virtual TRefArray* GetClusters() = 0;
+
+ // Return true if cluster has energy > cut
+ virtual Bool_t CutMinEnergy(const AliESDCaloCluster & /*cluster*/) const { return true; }
+
+ // Cut on distance to bad channel
+ virtual Bool_t CutDistanceToBadChannel(const AliESDCaloCluster & /*cluster*/) const { return true; }
+
+ // Cut on track matching
+ virtual Bool_t CutTrackMatching(const AliESDCaloCluster& /*cluster*/, Double_t &/*r*/) const { return true; }
+
+protected:
+
+ const AliVEvent *fEvent; // Pointer to current event
+
+ AliAnalysisEtCuts *fCuts; // Pointer to the cuts object
+
+ Int_t fRunNumber;
+
+private:
+
+ AliAnalysisEtSelector(); // Prohibited
+ AliAnalysisEtSelector(const AliAnalysisEtSelector& other);// Prohibited
+ AliAnalysisEtSelector& operator=(const AliAnalysisEtSelector& other);// Prohibited
+ bool operator==(const AliAnalysisEtSelector& other) const;// Prohibited
+
+ ClassDef(AliAnalysisEtSelector, 1);
+};
+
+#endif // ALIANALYSISETSELECTOR_H
--- /dev/null
+#include "AliAnalysisEtSelectorPhos.h"
+#include "AliAnalysisEtCuts.h"
+#include "AliESDCaloCluster.h"
+#include <AliVEvent.h>
+#include "TRefArray.h"
+#include "AliPHOSGeometry.h"
+#include "TH2I.h"
+#include "TFile.h"
+#include "TMath.h"
+#include "AliLog.h"
+#include <iostream>
+ClassImp(AliAnalysisEtSelectorPhos)
+AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
+,fGeoUtils(0)
+,fBadMapM2(0)
+,fBadMapM3(0)
+,fBadMapM4(0)
+{
+
+}
+
+AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
+{
+
+}
+
+TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
+{
+ return 0;
+}
+
+Int_t AliAnalysisEtSelectorPhos::Init(Int_t runNumber)
+{
+ AliAnalysisEtSelector::Init(runNumber);
+
+ int res = LoadGeometry();
+ if(res) return -1;
+ return LoadBadMaps();
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster) const
+{
+ return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const
+{
+ Float_t gPos[3];
+ cluster.GetPosition(gPos);
+ Int_t relId[4];
+ TVector3 glVec(gPos);
+ fGeoUtils->GlobalPos2RelId(glVec, relId);
+
+ TVector3 locVec;
+ fGeoUtils->Global2Local(locVec, glVec, relId[0]);
+// std::cout << fGeoUtils << std::endl;
+ //std::cout << relId[0] << " " << cluster.IsPHOS() << std::endl;
+ //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
+ for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
+ {
+ for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
+ {
+ if (relId[0] == 3)
+ {
+ if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 3;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kFALSE;
+ }
+ }
+ }
+ if (relId[0] == 2)
+ {
+ if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 2;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kFALSE;
+ }
+ }
+ }
+ if (relId[0] == 1)
+ {
+ if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
+ {
+ Int_t tmpRel[4];
+ tmpRel[0] = 1;
+ tmpRel[1] = 0;
+ tmpRel[2] = x+1;
+ tmpRel[3] = z+1;
+
+ Float_t tmpX;
+ Float_t tmpZ;
+ fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
+
+ Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
+
+// Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
+ //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
+ if (distance < fCuts->GetPhosBadDistanceCut())
+ {
+// std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() << std::endl;
+ return kFALSE;
+ }
+ }
+ }
+
+ }
+ }
+
+ return kTRUE;
+
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& cluster, Double_t &r) const
+{
+
+
+ // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
+
+ Int_t nTracksMatched = cluster.GetNTracksMatched();
+ if(nTracksMatched == 0) return kFALSE;
+
+ Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
+ if(trackMatchedIndex < 0) return kTRUE;
+
+ AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
+ Double_t dx = cluster.GetTrackDx();
+ Double_t dz = cluster.GetTrackDz();
+ Double_t pt = track->Pt();
+ Int_t charge = track->Charge();
+
+ Double_t meanX=0;
+ Double_t meanZ=0.;
+ Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
+ 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
+ Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
+
+ Double_t mf = fEvent->GetMagneticField(); //Positive for ++ and negative for --
+
+ if(mf<0.){ //field --
+ meanZ = -0.468318 ;
+ if(charge>0)
+ meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
+ else
+ meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
+ }
+ else{ //Field ++
+ meanZ= -0.468318;
+ if(charge>0)
+ meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
+ else
+ meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
+ }
+
+ Double_t rz=(dz-meanZ)/sz ;
+ Double_t rx=(dx-meanX)/sx ;
+ r = TMath::Sqrt(rx*rx+rz*rz);
+ if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
+
+ return kTRUE;
+
+}
+
+int AliAnalysisEtSelectorPhos::LoadGeometry()
+{
+
+ fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
+ // ifstream f("badchannels.txt", ios::in);
+ return 0;
+}
+
+int AliAnalysisEtSelectorPhos::LoadBadMaps()
+{
+TFile *f = TFile::Open("badchannels.root", "READ");
+
+ if(!f)
+ {
+ std::cout << "Could not open badchannels.root" << std::endl;
+ return -1;
+ }
+
+ fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
+ if(fBadMapM2)
+ {
+ std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
+ }
+ fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
+ if(fBadMapM2)
+ {
+ std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
+ }
+
+ fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
+ if(fBadMapM4)
+ {
+ std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
+ }
+
+ return 0;
+
+
+
+}
--- /dev/null
+#ifndef ALIANALYSISETSELECTORPHOS_H
+#define ALIANALYSISETSELECTORPHOS_H
+
+#include "AliAnalysisEtSelector.h"
+
+class TH2I;
+class AliPHOSGeometry;
+
+class AliAnalysisEtSelectorPhos : public AliAnalysisEtSelector
+{
+
+public:
+
+ AliAnalysisEtSelectorPhos(AliAnalysisEtCuts *cuts);
+ virtual ~AliAnalysisEtSelectorPhos();
+
+ virtual TRefArray* GetClusters();
+ virtual Bool_t CutMinEnergy(const AliESDCaloCluster& cluster) const;
+ virtual Bool_t CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const;
+ virtual Bool_t CutTrackMatching(const AliESDCaloCluster& cluster, Double_t &r) const;
+
+ virtual Int_t Init(int runNumber);
+
+private:
+
+
+ int LoadGeometry();
+ int LoadBadMaps();
+
+ AliPHOSGeometry *fGeoUtils;
+
+ TH2I *fBadMapM2; // Bad map
+ TH2I *fBadMapM3; // Bad map
+ TH2I *fBadMapM4; // Bad map
+
+ AliAnalysisEtSelectorPhos();
+ AliAnalysisEtSelectorPhos(const AliAnalysisEtSelectorPhos& other); // Prohibited
+ AliAnalysisEtSelectorPhos& operator=(const AliAnalysisEtSelectorPhos& other); // Prohibited
+ bool operator==(const AliAnalysisEtSelectorPhos& other) const; // Prohibited
+
+ ClassDef(AliAnalysisEtSelectorPhos, 1);
+};
+
+#endif // ALIANALYSISETSELECTORPHOS_H
,fEsdtrackCutsTPC(0)
,fEsdtrackCutsITS(0)
,fOutputList(0)
+ ,fPhysSelTaskName("physSelTask")
,fCentSelTaskName("centralityTask")
,fIsMc(isMc)
,fCurrentRunNum(-1)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff