}
// calculate ET
- Double_t etDep = CalculateTransverseEnergy(caloCluster);
+ Double_t etDep = CalculateTransverseEnergy(*caloCluster);
// All clusters
//fHistAllRecEtaEDepETDep->Fill(caloE,caloPos.Eta(),etDep);
#include "TString.h"
#include "AliCentrality.h"
#include "AliAnalysisEtSelector.h"
-
- //#include "THnSparse.h"
+#include "AliAnalysisEtTrackMatchCorrections.h"
+#include "AliAnalysisEtRecEffCorrection.h"
+#include "TFile.h"
+#include "TVector3.h"
using namespace std;
ClassImp(AliAnalysisEt);
AliAnalysisEt::AliAnalysisEt() : AliAnalysisEtCommon()
+ ,fTmCorrections(0)
+ ,fReCorrections(0)
,fEventSummaryTree(0)
,fAcceptedTree(0)
,fDepositTree(0)
void AliAnalysisEt::Init()
{// clear variables, set up cuts and PDG info
AliAnalysisEtCommon::Init();
+ if(ReadCorrections("calocorrections.root") != 0)
+ {
+ // Shouldn't do this, why oh why are exceptions not allowed?
+ exit(-1);
+ }
ResetEventValues();
}
return;
}
-Double_t AliAnalysisEt::CalculateTransverseEnergy(AliESDCaloCluster* cluster)
-{ // based on cluster energy and cluster pos
-
- Float_t pos[3];
- cluster->GetPosition(pos);
- TVector3 cp(pos);
- Double_t corrEnergy = 0;
-
- if(cluster->E() < 1.5)
+Int_t AliAnalysisEt::ReadCorrections(TString filename)
+{
+ TFile *f = TFile::Open(filename, "READ");
+ if(f)
{
- corrEnergy =cluster->E()/(0.51 + 0.02*cluster->E());
- }
- else
- {
- corrEnergy =cluster->E()/(0.51 + 0.02*1.5);
+ TString det = "Phos";
+ if(fHistogramNameSuffix.Contains("Emcal"))
+ {
+ det = "Emcal";
+ }
+ TString name = "TmCorrections" + det;
+ std::cout << name << std::endl;
+ fTmCorrections = dynamic_cast<AliAnalysisEtTrackMatchCorrections*>(f->Get(name));
+ if(!fTmCorrections)
+ {
+ Printf("Could not load TM corrections");
+ return -1;
+ }
+ name = "ReCorrections" + det;
+ fReCorrections = dynamic_cast<AliAnalysisEtRecEffCorrection*>(f->Get(name));
+ if(!fReCorrections)
+ {
+ Printf("Could not load rec eff corrections");
+ return -1;
+ }
+ return 0;
}
- //std::cout << "Original energy: " << cluster->E() << ", corrected energy: " << corrEnergy << std::endl;
+ return -1;
+}
+
+Double_t AliAnalysisEt::CalculateTransverseEnergy(const AliESDCaloCluster& cluster)
+{
+ Float_t pos[3];
+ cluster.GetPosition(pos);
+ TVector3 cp(pos);
+ Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E());
- return corrEnergy * TMath::Sin(cp.Theta());
+// std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl;
+ return TMath::Sin(cp.Theta())*corrEnergy;
}
#include "THnSparse.h"
#include "AliESDCaloCluster.h"
#include "AliAnalysisEtCuts.h"
+#include "AliAnalysisEtTrackMatchCorrections.h"
#include <vector>
#include "Rtypes.h"
+class AliAnalysisEtRecEffCorrection;
+class AliAnalysisEtTrackMatchCorrections;
class AliAnalysisEtSelector;
class AliCentrality;
class TString;
}
/** Get contribution from non-removed charged particles */
- virtual Double_t GetChargedContribution(Int_t /*clusterMultiplicity*/) {
- return 0;
+ Double_t GetChargedContribution(Int_t clusterMultiplicity) {
+ return fTmCorrections->ChargedContr(clusterMultiplicity);
}
/** Get contribution from non-removed neutral particles */
- virtual Double_t GetNeutralContribution(Int_t /*clusterMultiplicity*/) {
- return 0;
+ Double_t GetNeutralContribution(Int_t clusterMultiplicity) {
+ return fTmCorrections->NeutralContr(clusterMultiplicity);
}
/** Get contribution from removed gammas */
- virtual Double_t GetGammaContribution(Int_t /*clusterMultiplicity*/) {
- return 0;
+ Double_t GetGammaContribution(Int_t clusterMultiplicity) {
+ return fTmCorrections->GammaContr(clusterMultiplicity);
}
/** Get contribution from secondaries */
- virtual Double_t GetSecondaryContribution(Int_t /*clusterMultiplicity*/) {
- return 0;
+ Double_t GetSecondaryContribution(Int_t clusterMultiplicity) {
+ return fTmCorrections->SecondaryContr(clusterMultiplicity);
}
void MakeSparseHistograms() {
}
AliAnalysisEtCuts * GetCuts() const { return fCuts; }
+
+
+ // Read in corrections
+ Int_t ReadCorrections(TString filename); // Read in corrections
+
protected:
//AliAnalysisEtCuts *fCuts; // keeper of basic cuts
- Double_t CalculateTransverseEnergy(AliESDCaloCluster *cluster);
+
+ // Return corrected cluster E_T
+ Double_t CalculateTransverseEnergy(const AliESDCaloCluster &cluster);
+
+ // Track matching (hadrdonic contamination) corrections
+ AliAnalysisEtTrackMatchCorrections *fTmCorrections;
+
+ // Reconstruction efficiency corrections
+ AliAnalysisEtRecEffCorrection *fReCorrections;
TTree *fEventSummaryTree; // Contains event level information
AliAnalysisEtSelector *fSelector; // Selector class
private:
+
+
+
//Declare private to avoid compilation warning
AliAnalysisEt & operator = (const AliAnalysisEt & g) ;//cpy assignment
AliAnalysisEt(const AliAnalysisEt & g) ; // cpy ctor
Int_t AliAnalysisEtCommon::fgAntiProtonCode = -2212;
Int_t AliAnalysisEtCommon::fgLambdaCode = 3122;
Int_t AliAnalysisEtCommon::fgAntiLambdaCode = -3122;
+Int_t AliAnalysisEtCommon::fgK0Code = 311;
Int_t AliAnalysisEtCommon::fgK0SCode = 310;
Int_t AliAnalysisEtCommon::fgOmegaCode = 3334;
Int_t AliAnalysisEtCommon::fgAntiOmegaCode = -3334;
static Int_t fgAntiProtonCode;//pdg antiproton code
static Int_t fgLambdaCode;// pdg lambda code
static Int_t fgAntiLambdaCode;//pdg antilambda code
- static Int_t fgK0SCode;//pdg k0 short code
static Int_t fgOmegaCode;//pdg omega code
+ static Int_t fgK0SCode;//pdg k0 short code
+ static Int_t fgK0Code;//pdg k0 code
static Int_t fgAntiOmegaCode;//pdg anti-omega code
static Int_t fgXi0Code;//pdg xi-0 code
static Int_t fgAntiXi0Code;//pdg anti-xi0 code
,fPhosTrackDistanceCut(10.0)
,fPhosTrackDxCut(8.0)
,fPhosTrackDzCut(3.0)
- ,fPhosTrackRCut(2.0)
+ ,fPhosTrackRCut(5.0)
,fPhosBadDistanceCut(3.0)
,fGeometryPhosEtaAccCut(0.12)
,fReconstructedPidCut(0.0)
//
,fReconstructedPhosClusterType(-1)
- ,fReconstructedPhosClusterEnergyCut(0.15)
+ ,fReconstructedPhosClusterEnergyCut(0.25)
,fReconstructedPhosSingleCellEnergyCut(0.5)
,fReconstructedPhosTrackDistanceTightCut(3.0)
,fReconstructedPhosTrackDistanceMediumCut(5.0)
,fHistNbinsParticlePt(200)
,fHistMinParticlePt(0)
,fHistMaxParticlePt(20)
+
+ ,fPrimaryVertexCutXY(4.0)
+ ,fPrimaryVertexCutZ(20.0)
{ // ctor
}
//_________________________________________________________________________
#include "TNamed.h"
+#include <iostream>
class AliAnalysisEtCuts : public TNamed
{
Short_t GetDetectorPhos() const { return fgkDetectorPhos; }
Short_t GetDetectorEmcal() const { return fgkDetectorEmcal; }
+
+ Double_t GetPrimaryVertexCutXY() const { return fPrimaryVertexCutXY; }
+ Double_t GetPrimaryVertexCutZ() const { return fPrimaryVertexCutZ; }
+
// Setters
// Common
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 SetPhosTrackRCut(const Double_t val) { std::cout << "Setting: " << val << std::endl; fPhosTrackRCut = val; }
void SetPhosBadDistanceCut(const Double_t val) { fPhosBadDistanceCut = val; }
void SetHistMinParticlePt(const Double_t val) { fHistMinParticlePt = val; }
void SetHistMaxParticlePt(const Double_t val) { fHistMaxParticlePt = val; }
-
+ void SetPrimaryVertexCutXY(const Double_t val) { fPrimaryVertexCutXY = val; }
+ void SetPrimaryVertexCutZ(const Double_t val) { fPrimaryVertexCutZ = val; }
+
+
+
protected:
// Common
// Detector definition
static const Short_t fgkDetectorPhos = -1; // PHOS
static const Short_t fgkDetectorEmcal = 1; // EMCAL
+
+ Double_t fPrimaryVertexCutXY; // Cut to decide if particle is from primary vertex
+ Double_t fPrimaryVertexCutZ; // Cut to decide if particle is from primary vertex
+
private:
//Declare private to avoid compilation warning
,fPrimaryVy(0)
,fPrimaryVz(0)
,fPrimaryAccepted(0)
+ ,fPrimaryMatched(0)
,fDepositedCode(0)
+ ,fDepositedE(0)
,fDepositedEt(0)
,fDepositedCharge(0)
,fDepositedVx(0)
,fDepositedVy(0)
,fDepositedVz(0)
+ ,fSecondary(kFALSE)
+ ,fReconstructedE(0)
+ ,fReconstructedEt(0)
+ ,fTotPx(0)
+ ,fTotPy(0)
+ ,fTotPz(0)
+ ,fClusterMult(0)
,fHistDecayVertexNonRemovedCharged(0)
,fHistDecayVertexRemovedCharged(0)
,fHistDecayVertexNonRemovedNeutral(0)
,fHistPiPlusMultAcc(0)
,fHistPiMinusMultAcc(0)
,fHistPiZeroMultAcc(0)
- ,fPiPlusMult(0)
- ,fPiMinusMult(0)
+// ,fPiPlusMult(0)
+// ,fPiMinusMult(0)
,fPiZeroMult(0)
,fPiPlusMultAcc(0)
,fPiMinusMultAcc(0)
,fChargedRemoved(0)
,fChargedNotRemoved(0)
,fNeutralNotRemoved(0)
+ ,fGammaRemoved(0)
+ ,fSecondaryNotRemoved(0)
,fEnergyNeutralRemoved(0)
,fEnergyChargedRemoved(0)
,fEnergyChargedNotRemoved(0)
,fEnergyNeutralNotRemoved(0)
+ ,fEnergyGammaRemoved(0)
,fNClusters(0)
,fTotNeutralEtAfterMinEnergyCut(0)
-
+ ,fHistGammasFound(0)
+ ,fHistGammasGenerated(0)
{
}
{ // analyse MC event
ResetEventValues();
+
fPiPlusMult = 0;
fPiMinusMult = 0;
fPiZeroMult = 0;
{
TParticle *part = stack->Particle(iPart);
+
+
if (!part)
{
Printf("ERROR: Could not get particle %d", iPart);
continue;
}
-
TParticlePDG *pdg = part->GetPDG(0);
if (!pdg)
Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
Int_t code = pdg->PdgCode();
-
+
+ if(stack->IsPhysicalPrimary(iPart))
+ {
+ fTotPx += part->Px();
+ fTotPy += part->Py();
+ fTotPz += part->Pz();
+ }
+
// Check for reasonable (for now neutral and singly charged) charge on the particle
- //TODO:Maybe not only singly charged?
if (fSelector->CutNeutralMcParticle(iPart,*stack,*pdg))
{
fMultiplicity++;
- //PrintFamilyTree(iPart, stack);
+// PrintFamilyTree(iPart, stack);
//
if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
{
if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
{
-
+ // PrintFamilyTree(iPart,stack);
//Printf("Gamma, phi: %f, eta: %f, phi cut min: %f, phi cut max: %f, eta cut: %f", part->Phi(), part->Eta(), fPhiMinCutAcc, fPhiMaxCutAcc, fEtaCutAcc);
//if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
//if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
- if(fSelector->CutGeometricalAcceptance(*part))
+ if(fSelector->CutGeometricalAcceptance(*part) )
{
fNeutralMultiplicity++;
fTotNeutralEt += et;
}
}
-
+ // std::cout << "Total: p_x = " << fTotPx << ", p_y = " << fTotPy << ", p_z = " << fTotPz << std::endl;
fTotEt = fTotChargedEt + fTotNeutralEt;
//fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
//std::cout << "Event done! # of particles: " << partCount << std::endl;
AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
if (!mcEvent || !realEvent) {
AliFatal("ERROR: mcEvent or realEvent does not exist");
- return 0;
+
}
+ std::vector<Int_t> foundGammas;
+
fSelector->SetEvent(realEvent);
AnalyseEvent(ev);
TRefArray *caloClusters = fSelector->GetClusters();
Int_t nCluster = caloClusters->GetEntries();
+ fClusterMult = nCluster;
fNClusters = 0;
// loop the clusters
for (int iCluster = 0; iCluster < nCluster; iCluster++ )
Int_t cf = 0;
AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
//Float_t caloE = caloCluster->E();
- fCutFlow->Fill(cf++);
- if(!fSelector->CutDistanceToBadChannel(*caloCluster)) continue;
fNClusters++;
- UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
+ const UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
+ const UInt_t childPart = iPart;
TParticle *part = stack->Particle(iPart);
if (!part)
if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
{
primIdx = GetPrimMother(iPart, stack);
+ if(primIdx != stack->GetPrimary(iPart))
+ {
+ //std::cout << primIdx << " = " << stack->GetPrimary(iPart) << std::endl;
+ //PrintFamilyTree(iPart, stack);
+ }
+
if(primIdx < 0)
{
std::cout << "What!? No primary?" << std::endl;
}
} // end of primary particle check
+
const int primCode = stack->Particle(primIdx)->GetPdgCode();
TParticlePDG *pdg = part->GetPDG();
if (!pdg)
}
Int_t code = pdg->PdgCode();
- Double_t clEt = CalculateTransverseEnergy(caloCluster);
+ if(primCode == fgGammaCode)
+ {
+
+
+ if(fSelector->CutDistanceToBadChannel(*caloCluster))//&&fSelector->CutGeometricalAcceptance(*(stack->Particle(primIdx))))
+ {
+// std::cout << "Gamma primary: " << primIdx << std::endl;
+ foundGammas.push_back(primIdx);
+ }
+
+ }
+ fCutFlow->Fill(cf++);
+ if(!fSelector->CutDistanceToBadChannel(*caloCluster)) continue;
+ Double_t clEt = CalculateTransverseEnergy(*caloCluster);
// if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
if(!fSelector->CutMinEnergy(*caloCluster)) continue;
fCutFlow->Fill(cf++);
Float_t pos[3];
-
+ //PrintFamilyTree(
caloCluster->GetPosition(pos);
TVector3 cp(pos);
fPrimaryVz = primPart->Vz();
fPrimaryAccepted = false;
+ fPrimaryMatched = false;
fDepositedCode = part->GetPdgCode();
- fDepositedEt = caloCluster->E() * TMath::Sin(cp.Eta());
+ fDepositedE = part->Energy();
+ fDepositedEt = part->Energy()*TMath::Sin(part->Theta());
fDepositedCharge = part->GetPDG()->Charge();
fDepositedVx = part->Vx();
fDepositedVy = part->Vy();
fDepositedVz = part->Vz();
+
+ //fSecondary = fSelector->FromSecondaryInteraction(*primPart, *stack);
+ fSecondary =fSelector->FromSecondaryInteraction(*part, *stack);
+ if(fSecondary)
+ {
+ //std::cout << "Have secondary!" << std::endl;
+ //PrintFamilyTree(iPart, stack);
+ }
+ fReconstructedE = caloCluster->E();
+ fReconstructedEt = caloCluster->E()*TMath::Sin(cp.Theta());
//if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
- pdg = part->GetPDG(0);
- //if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
+
+ pdg = primPart->GetPDG(0);
+ code = primPart->GetPdgCode();
+ //if (TMath::Abs(caloCluster->GetTrackDx()) < 5 && TMath::Abs(caloCluster->GetTrackDz()) < 5)
if(!fSelector->CutTrackMatching(*caloCluster))
- //if(caloCluster->GetTrackMatchedIndex() > -1 && (fCuts->GetPhosTrackRCut() > TestCPV(caloCluster->GetTrackDx(), caloCluster->GetTrackDz(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Pt(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Charge(), ev)))
{
-
+ fPrimaryMatched = true;
if (pdg->Charge() != 0)
{
- //std::cout << "Removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
fChargedRemoved++;
- //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());
fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
- if(code == fgProtonCode)
- {
- fProtonRemovedEt += clEt;
- fProtonRemovedMult++;
-
- }
- else if(code == fgAntiProtonCode)
- {
- fAntiProtonRemovedEt += clEt;
- fAntiProtonRemovedMult++;
- }
- else if(code == fgPiPlusCode)
- {
- fPiPlusRemovedEt += clEt;
- fPiPlusRemovedMult++;
- }
- else if(code == fgPiMinusCode)
- {
- fPiMinusRemovedEt += clEt;
- fPiMinusRemovedMult++;
- }
- else if(code == fgKPlusCode)
- {
- fKPlusRemovedEt += clEt;
- fKPlusRemovedMult++;
- }
- else if(code == fgKMinusCode)
- {
- fKMinusRemovedEt += clEt;
- fKMinusRemovedMult++;
- }
- else if(code == fgMuMinusCode)
- {
- fMuMinusRemovedEt += clEt;
- fMuMinusRemovedMult++;
- }
- else if(code == fgMuPlusCode)
- {
- fMuPlusRemovedEt += clEt;
- fMuPlusRemovedMult++;
- }
- else if(code == fgEMinusCode)
- {
- fEMinusRemovedEt += clEt;
- fEMinusRemovedMult++;
- }
- else if(code == fgEPlusCode)
- {
- fEPlusRemovedEt += clEt;
- fEPlusRemovedMult++;
- }
- else
- {
- Printf("Removed: %d, with E_T: %f", code, clEt);
- }
- if (primCode == fgGammaCode)
- {
- fGammaRemovedEt += clEt;
- fGammaRemovedMult++;
- }
- else if (primCode == fgPi0Code)
- {
- fPi0RemovedEt += clEt;
- fPi0RemovedMult++;
- }
- else if (primCode == fgEtaCode)
- {
- fPi0RemovedEt += clEt;
- fPi0RemovedMult++;
- }
- else if(primCode == fgK0LCode)
- {
- fK0lRemovedEt += clEt;
- fK0lRemovedMult++;
- }
- else if(primCode == fgK0SCode)
- {
- fK0sRemovedEt += clEt;
- fK0sRemovedMult++;
- }
-
+
}
else
{
- //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 += 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 (code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
- {
- fEtRemovedGammas += clEt;
- fMultRemovedGammas++;
- }
- else if (code == fgNeutronCode)
- {
- fEtRemovedNeutrons += clEt;
- fMultRemovedNeutrons++;
- fNeutronRemovedEt += clEt;
- fNeutronRemovedMult++;
- }
- else if (code == fgAntiNeutronCode)
- {
- fEtRemovedAntiNeutrons += clEt;
- fMultRemovedAntiNeutrons++;
- fAntiNeutronRemovedEt += clEt;
- fAntiNeutronRemovedMult++;
- }
- if (primCode == fgGammaCode)
- {
- fGammaRemovedEt += clEt;
- fGammaRemovedMult++;
- }
- else if(primCode == fgPi0Code)
- {
- fPi0RemovedEt += clEt;
- fPi0RemovedMult++;
- }
- else if(primCode == fgEtaCode)
- {
- fPi0RemovedEt += clEt;
- fPi0RemovedMult++;
- }
- else if(primCode == fgK0LCode)
- {
- fK0lRemovedEt += clEt;
- fK0lRemovedMult++;
- }
- else if(primCode == fgK0SCode)
- {
- fK0sRemovedEt += clEt;
- fK0sRemovedMult++;
- }
+ if(code==fgGammaCode)
+ {
+ fGammaRemoved++;
+ fGammaRemovedEt+=clEt;
+
+ }
+ else
+ {
+ fNeutralRemoved++;
+ fEnergyNeutralRemoved += clEt;
+ fHistRemovedOrNot->Fill(1.0, fCentClass);
+ fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
+ }
}
}
else
{
+ fPrimaryAccepted = true;
+ fPrimaryMatched = false;
fCutFlow->Fill(cf++);
if (pdg->Charge() != 0)
{
- //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 += 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 (code == fgProtonCode)
- {
- //std::cout << clEt << std::endl;
- fEtNonRemovedProtons += clEt;
- fMultNonRemovedProtons++;
- fProtonEt += clEt;
- fProtonMult++;
- }
- else if (code == fgAntiProtonCode)
- {
- //std::cout << clEt << std::endl;
- fEtNonRemovedAntiProtons += clEt;
- fMultNonRemovedAntiProtons++;
- fAntiProtonEt += clEt;
- fAntiProtonMult++;
- }
- else if (code == fgPiPlusCode)
- {
- //std::cout << "PI+" << clEt << std::endl;
- fEtNonRemovedPiPlus += clEt;
- fMultNonRemovedPiPlus++;
- fPiPlusEt += clEt;
- fPiPlusMult++;
- }
- else if (code == fgPiMinusCode)
- {
- // std::cout << "PI-" << clEt << std::endl;
- fEtNonRemovedPiMinus += clEt;
- fMultNonRemovedPiMinus++;
- fPiMinusEt += clEt;
- fPiMinusMult++;
- }
- else if (code == fgKPlusCode)
- {
- //std::cout << clEt << std::endl;
- fEtNonRemovedKaonPlus += clEt;
- fMultNonRemovedKaonPlus++;
- fKPlusEt += clEt;
- fKPlusMult++;
- }
- else if (code == fgKMinusCode)
- {
- //std::cout << clEt << std::endl;
- fEtNonRemovedKaonMinus += clEt;
- fMultNonRemovedKaonMinus++;
- fKMinusEt += clEt;
- fKMinusMult++;
- }
- else if (code == fgEPlusCode)
- {
- //std::cout << clEt << std::endl;
- if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 440)
- {
- fEtNonRemovedPositrons += clEt;
- fMultNonRemovedPositrons++;
- fEPlusEt += clEt;
- fEPlusMult++;
- }
- }
- else if (code == fgEMinusCode)
- {
- //std::cout << clEt << std::endl;
- if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 440)
- {
- fEtNonRemovedElectrons += clEt;
- fMultNonRemovedElectrons++;
- fEMinusEt += clEt;
- fEMinusMult++;
- }
- }
- else if (code == fgMuPlusCode)
- {
- fEtNonRemovedMuPlus += clEt;
- fMultNonRemovedMuPlus++;
- fMuPlusEt += clEt;
- fMuPlusMult++;
- }
- else if (code == fgMuMinusCode)
- {
- fEtNonRemovedMuMinus += clEt;
- fMultNonRemovedMuMinus++;
- fMuMinusEt += clEt;
- fMuMinusMult++;
- }
- if (primCode == fgGammaCode)
- {
- fGammaEt += clEt;
- fGammaMult++;
- }
- else if(primCode == fgPi0Code)
- {
- fPi0Et += clEt;
- fPi0Mult++;
- }
- else if(primCode == fgEtaCode)
- {
- fPi0Et += clEt;
- fPi0Mult++;
- }
- else if(primCode == fgK0LCode)
- {
- fK0lEt += clEt;
- fK0lMult++;
- }
- else if(primCode == fgK0SCode)
- {
- fK0sEt += clEt;
- fK0sMult++;
- }
+ if(fSecondary)
+ {
+ fSecondaryNotRemoved++;
+ }
+ else
+ {
+ fChargedNotRemoved++;
+ fEnergyChargedNotRemoved += clEt;
+ fHistRemovedOrNot->Fill(2.0, fCentClass);
+ fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
+ PrintFamilyTree(iPart, stack);
+ }
}
else
{
- fPrimaryAccepted = true;
- //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 += clEt;
- fHistRemovedOrNot->Fill(3.0, fCentClass);
- fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
-
- if(TMath::Abs(part->Vx()) < 1.0 && TMath::Abs(part->Vy()) < 1.0 && TMath::Abs(part->Vz()) < 20 && fSelector->IsEmEtParticle(primCode))
+ if(!fSecondary)
+ {
+ fNeutralNotRemoved++;
+ fEnergyNeutralNotRemoved += clEt;
+ fHistRemovedOrNot->Fill(3.0, fCentClass);
+ fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
+ }
+ else
+ {
+// PrintFamilyTree(iPart, stack);
+ }
+ // if(TMath::Abs(part->Vx()) < 1.0 && TMath::Abs(part->Vy()) < 1.0 && TMath::Abs(part->Vz()) < 20 && fSelector->IsEmEtParticle(primCode))
+ if(fSecondary && fSelector->IsEmEtParticle(primCode))
{
- fTotEtWithSecondaryRemoved += clEt;
+ fTotEtSecondaryFromEmEtPrimary += clEt;
+ fSecondaryNotRemoved++;
}
else if(fSelector->IsEmEtParticle(primCode))
{
- fTotEtSecondaryFromEmEtPrimary += clEt;
+ fTotEtWithSecondaryRemoved += clEt;
+// PrintFamilyTree(iPart, stack);
+
}
- else {
- fTotEtSecondary += clEt;
+ else
+ {
+ //fTotEtSecondary += clEt;
+// PrintFamilyTree(iPart,stack);
}
//code = stack->Particle(primIdx)->GetPdgCode();
- if (code == fgGammaCode)
+ if (code == fgGammaCode && !fSecondary)
{
fEtNonRemovedGammas += clEt;
fMultNonRemovedGammas++;
+ fNeutralNotRemoved--;
+ fEnergyNeutralNotRemoved -= clEt;
+ //PrintFamilyTree(iPart, stack);
// if (pdgMom)
// {
// if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
// }
// }
}
- 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 (code == fgAntiNeutronCode)
- {
- fEtNonRemovedAntiNeutrons += clEt;
- fMultNonRemovedAntiNeutrons++;
- }
- //else if (code == fgK0LCode || pdgMom->PdgCode() == fgK0SCode)
- else if (code == fgK0SCode)
- {
- //std::cout << "K0 with energy: " << clEt << std::endl;
- fEtNonRemovedK0S += clEt;
- fMultNonRemovedK0s++;
- }
- 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): " << code << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
- if (primCode == fgGammaCode)
- {
- fGammaEt += clEt;
- fGammaMult++;
- }
- else if(primCode == fgPi0Code)
- {
- fPi0Et += clEt;
- fPi0Mult++;
- }
- else if(primCode == fgEtaCode)
- {
- fPi0Et += clEt;
- fPi0Mult++;
- }
- else if(primCode == fgK0LCode)
- {
- fK0lEt += clEt;
- fK0lMult++;
- }
- else if(primCode == fgK0SCode)
- {
- fK0sEt += clEt;
- fK0sMult++;
- }
}
}
fPrimaryTree->Fill();
} // end of loop over clusters
+
+ std::sort(foundGammas.begin(), foundGammas.end());
+ for (Int_t iPart = 0; iPart < stack->GetNtrack(); iPart++)
+ {
+
+ if(!stack->IsPhysicalPrimary(iPart)) continue;
+
+ TParticle *part = stack->Particle(iPart);
+
+ if (!part)
+ {
+ Printf("ERROR: Could not get particle %d", iPart);
+ continue;
+ }
+ TParticlePDG *pdg = part->GetPDG(0);
+
+ if (!pdg)
+ {
+ Printf("ERROR: Could not get particle PDG %d", iPart);
+ continue;
+ }
+
+ if(pdg->PdgCode()==fgGammaCode && fSelector->CutGeometricalAcceptance(*part))// TMath::Abs(part->Eta()) < 0.12)
+ {
+ fHistGammasGenerated->Fill(part->Energy());
+// std::cout << "Searching for gammas..." << std::endl;
+ if(std::binary_search(foundGammas.begin(),foundGammas.end(),iPart))
+ {
+ fHistGammasFound->Fill(part->Energy());
+ }
+// for(int i = 0; i < foundGammas.size(); i++)
+// {
+// // std::cout << iPart << std::endl;
+// if(foundGammas[i] == iPart)
+// {
+// fHistGammasFound->Fill(part->Energy());
+// std::cout << "Match!" << std::endl;
+// break;
+// }
+// }
+ }
+
+ }
//std::cout << "Distance cut: " << fTrackDistanceCut << std::endl;
//std::cout << "Number of removed neutrals: " << fNeutralRemoved << std::endl;
//std::cout << "Number of removed charged: " << fChargedRemoved << std::endl;
fChargedRemoved = 0;
fNeutralNotRemoved = 0;
fNeutralRemoved = 0;
-
+ fGammaRemoved = 0;
+ fSecondaryNotRemoved = 0;
fTrackMultInAcc = 0;
fTotNeutralEtAfterMinEnergyCut = 0;
-
+
+ fSecondaryNotRemoved = 0;
+
+ fTotPx = 0;
+ fTotPy = 0;
+ fTotPz = 0;
+
+
}
void AliAnalysisEtMonteCarlo::CreateHistograms()
fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
+ fEventSummaryTree->Branch("fSecondaryNotRemoved", &fSecondaryNotRemoved, "fSecondaryNotRemoved/I");
+ fEventSummaryTree->Branch("fChargedNotRemoved", &fChargedNotRemoved, "fChargedNotRemoved/I");
+ fEventSummaryTree->Branch("fNeutralNotRemoved", &fNeutralNotRemoved, "fNeutralNotRemoved/I");
+ fEventSummaryTree->Branch("fChargedRemoved", &fChargedRemoved, "fChargedRemoved/I");
+ fEventSummaryTree->Branch("fNeutralRemoved", &fNeutralRemoved, "fNeutralRemoved/I");
+ fEventSummaryTree->Branch("fGammaRemoved", &fGammaRemoved, "fGammaRemoved/I");
+ fEventSummaryTree->Branch("fTotPx", &fTotPx, "fTotPx/D");
+ fEventSummaryTree->Branch("fTotPy", &fTotPy, "fTotPy/D");
+ fEventSummaryTree->Branch("fTotPz", &fTotPz, "fTotPz/D");
+// fEventSummaryTree->Branch("f
}
//fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
+ fPrimaryTree->Branch("fPrimaryMatched", &fPrimaryMatched, "fPrimaryMatched/B");
fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
+ fPrimaryTree->Branch("fDepositedE", &fDepositedE, "fDepositedE/D");
fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
+ fPrimaryTree->Branch("fSecondary", &fSecondary, "fSecondary/I");
+
+
+ fPrimaryTree->Branch("fReconstructedE", &fReconstructedE, "fReconstructedE/D");
+ fPrimaryTree->Branch("fReconstructedEt", &fReconstructedEt, "fReconstructedEt/D");
+
+ fPrimaryTree->Branch("fClusterMult", &fClusterMult, "fClusterMult/I");
+
+
}
+ fHistGammasFound = new TH1F("fHistGammasFound", "fHistGammasFound",200, 0, 10);
+ fHistGammasGenerated = new TH1F("fHistGammasGenerated", "fHistGammasGenerated",200, 0, 10);
}
void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
list->Add(fHistPiPlusMultAcc);
list->Add(fHistPiMinusMultAcc);
list->Add(fHistPiZeroMultAcc);
+
+ list->Add(fHistGammasFound);
+ list->Add(fHistGammasGenerated);
}
+fMultNonRemovedProtons);
fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
- fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
+ fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
fMultRemovedGammas);
+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
- fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
+ fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas+fK0sMult);
fHistClusterMultvsRemovedGamma->Fill(fNClusters,
fMultRemovedGammas);
Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
{ // print family tree
TParticle *part = stack->Particle(partIdx);
- if(part->GetPdgCode() == fgK0SCode)
+// if(part->GetPdgCode() == fgK0SCode)
{
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;
return 0;
}
TParticle *mother = stack->Particle(mothIdx);
- if(mother->GetPdgCode() == fgK0SCode)
+// if(mother->GetPdgCode() == fgK0SCode)
{
- std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
+ //std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
+ std::cout << tabs << "Index: " << mothIdx << std::endl;
+ std::cout << tabs << "Primary: " << stack->IsPhysicalPrimary(mothIdx) << std::endl;
std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
+ if(mother->GetFirstMother() >= 0)
+ {
+ std::cout << tabs << "Mother(s): " << stack->Particle(mother->GetFirstMother())->GetPdgCode();
+ if(mother->GetSecondMother() >= 0) std::cout << ", " << stack->Particle(mother->GetSecondMother())->GetPdgCode();
+ std::cout << std::endl;
+ }
std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
}
if(mother->GetPdgCode() == fgK0SCode)
{ // get primary mother
if(partIdx >= 0)
{
+ //return stack->GetPrimary(partIdx);
+
Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
if(mothIdx < 0) return -1;
TParticle *mother = stack->Particle(mothIdx);
{ // get K0 in family
if(partIdx >= 0)
{
- if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
+ if(stack->Particle(partIdx)->GetPdgCode() == fgPi0Code) 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)
+ if(mother->GetPdgCode() == fgPi0Code)
{
return mothIdx;
}
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);
-
+
+
+
protected:
Double_t fImpactParameter; // b(fm), for Hijing; 0 otherwise
Int_t fNcoll; // Ncoll, for Hijing; 1 otherwise
Int_t fNpart; // Ncoll, for Hijing; 2 otherwise
-
+
TTree *fPrimaryTree; // Tree holding info on primaries
Double_t fTotEtWithSecondaryRemoved; // enter comment here
Double_t fPrimaryVz; // enter comment here
Bool_t fPrimaryAccepted; // enter comment here
- Int_t fDepositedCode; // enter comment here
- Double_t fDepositedEt; // enter comment here
+ Bool_t fPrimaryMatched;
+ Int_t fDepositedCode; // enter comment here Double_t fDepositedEt; // enter comment here
+ Double_t fDepositedE;
+ Double_t fDepositedEt;
Int_t fDepositedCharge; // enter comment here
Double_t fDepositedVx; // enter comment here
Double_t fDepositedVy; // enter comment here
Double_t fDepositedVz; // enter comment here
+
+ Bool_t fSecondary;
+
+ Double_t fReconstructedE;
+ Double_t fReconstructedEt;
+ Double_t fTotPx;
+ Double_t fTotPy;
+ Double_t fTotPz;
+
+ Int_t fClusterMult;
+
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
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
TH1F *fHistPiMinusMultAcc; // enter comment here
TH1F *fHistPiZeroMultAcc; // enter comment here
- // DS: pi+- mult already defined in base class..
- Int_t fPiPlusMult; // enter comment here
- Int_t fPiMinusMult; // enter comment here
+ // Int_t fPiPlusMult; // enter comment here
+ // Int_t fPiMinusMult; // enter comment here
Int_t fPiZeroMult; // enter comment here
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
+ Int_t fGammaRemoved; // number of gammas removed
+
+ Int_t fSecondaryNotRemoved;
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
+ Double_t fEnergyGammaRemoved; // energy of gammas that were removed
Int_t fNClusters; // Number of clusters in event
Double_t fTotNeutralEtAfterMinEnergyCut; // enter comment here
+ TH1F *fHistGammasFound;
+ TH1F *fHistGammasGenerated;
+
+
private:
//Declare it private to avoid compilation warning
--- /dev/null
+//_________________________________________________________________________
+// Utility Class for transverse energy studies
+// Selection class for EMCAL
+//
+//*-- Authors: Oystein Djuvsland (Bergen)
+//_________________________________________________________________________
+
+
+#include "AliAnalysisEtRecEffCorrection.h"
+
+ClassImp(AliAnalysisEtRecEffCorrection);
+
+AliAnalysisEtRecEffCorrection::AliAnalysisEtRecEffCorrection() : TNamed("RecEff","RecEff")
+ ,fEnergyCorrection("ReCorr","pol1",0.01)
+ ,fMaxEnergy(0)
+{}
+
+AliAnalysisEtRecEffCorrection::AliAnalysisEtRecEffCorrection(TString name, const TF1 &correction, const Double_t maxEnergy) : TNamed(name, name)
+ ,fEnergyCorrection(correction)
+ ,fMaxEnergy(maxEnergy)
+{}
+
+//! Copy constructor
+AliAnalysisEtRecEffCorrection::AliAnalysisEtRecEffCorrection(const AliAnalysisEtRecEffCorrection &obj) : TNamed(obj)
+ ,fEnergyCorrection(obj.fEnergyCorrection)
+ ,fMaxEnergy(obj.fMaxEnergy)
+{}
+
+//! Destructor
+AliAnalysisEtRecEffCorrection::~AliAnalysisEtRecEffCorrection()
+{}
+
+//! Assignment operator
+AliAnalysisEtRecEffCorrection& AliAnalysisEtRecEffCorrection::operator=(const AliAnalysisEtRecEffCorrection &other)
+{
+ if (this != &other)
+ {
+ fEnergyCorrection = other.fEnergyCorrection;
+ fMaxEnergy = other.fMaxEnergy;
+ }
+ return *this;
+}
+
+//! Equality operator
+bool AliAnalysisEtRecEffCorrection::operator==(const AliAnalysisEtRecEffCorrection &other) const
+{
+ if (this == &other) return true;
+ return false;
+// return (fMaxEnergy == other.fMaxEnergy && fEnergyCorrection == other.fEnergyCorrection);
+}
+
+Double_t AliAnalysisEtRecEffCorrection::CorrectedEnergy(Double_t energy)
+{
+
+ return fEnergyCorrection.Eval(energy);
+
+}
--- /dev/null
+//_________________________________________________________________________
+// Utility Class for transverse energy studies
+//
+//*-- Authors: Oystein Djuvsland (Bergen)
+//_________________________________________________________________________
+
+
+#ifndef ALIANALYSISETRECEFFCORRECTION_H
+#define ALIANALYSISETRECEFFCORRECTION_H
+
+#include "Rtypes.h"
+#include "TArrayD.h"
+#include "TNamed.h"
+#include "TF1.h"
+
+class AliAnalysisEtRecEffCorrection : public TNamed
+{
+
+public:
+
+//! Default constructor
+ AliAnalysisEtRecEffCorrection();
+
+//! Constructor
+ AliAnalysisEtRecEffCorrection(TString name, const TF1& correction, const Double_t maxEnergy);
+
+//! Copy constructor
+ AliAnalysisEtRecEffCorrection(const AliAnalysisEtRecEffCorrection &obj);
+
+//! Destructor
+ virtual ~AliAnalysisEtRecEffCorrection();
+
+//! Assignment operator
+ AliAnalysisEtRecEffCorrection& operator=(const AliAnalysisEtRecEffCorrection& other);
+
+//! Equality operator
+ bool operator==(const AliAnalysisEtRecEffCorrection &other) const;
+
+// Getters
+
+ TF1 EnergyCorrection() const {
+ return fEnergyCorrection;
+ }
+
+ Double_t MaxEnergy() const {
+ return fMaxEnergy;
+ }
+
+// Setters
+
+ void SetCorrections(const TF1 &corrections) {
+ fEnergyCorrection = corrections;
+ }
+
+ void SetMaxenergy(Double_t maxEnergy) {
+ fMaxEnergy = maxEnergy;
+ }
+
+
+ Double_t CorrectedEnergy(Double_t energy); // Calculate corrected cluster E_T
+
+private:
+
+ // Energy correction function
+ TF1 fEnergyCorrection;
+
+
+
+// MaxEnergy
+ Double_t fMaxEnergy; // MaxEnergy
+
+ ClassDef(AliAnalysisEtRecEffCorrection, 1);
+};
+
+#endif //ALIANALYSISETRECEFFCORRECTION_H
+
,fHistMuonEnergyDeposit(0)
,fHistRemovedEnergy(0)
,fGeomCorrection(1.0)
- ,fEMinCorrection(1.0/0.89)
+ ,fEMinCorrection(1.0/0.687)
,fRecEffCorrection(1.0)
,fClusterPosition(0)
,fHistChargedEnergyRemoved(0)
fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
- fTotNeutralEt += CalculateTransverseEnergy(cluster);
+ fTotNeutralEt += CalculateTransverseEnergy(*cluster);
fNeutralMultiplicity++;
}
fMultiplicity++;
fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
- Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
+ Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) + GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
fHistRemovedEnergy->Fill(removedEnergy);
fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
fTotEt = fTotChargedEt + fTotNeutralEt;
// Fill the histograms...0
FillHistograms();
- // std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
+ //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
return 0;
}
const Double_t kMEANGAMMA = 0.374;
*/
// Simulated case:
-//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);
+//corrEnergy =cluster->E()/(0.31 + 0.02*cluster->E());
+const Double_t kMEANCHARGED = 0.48/(0.31 + 0.00*0.48);
+const Double_t kMEANNEUTRAL = 0.53/(0.31 + 0.00*0.53);
+const Double_t kMEANGAMMA = 0.51/(0.31 + 0.00*0.31);
AliAnalysisEtReconstructedPhos::AliAnalysisEtReconstructedPhos() :
// fRemovedGammaContributionCorrectionParameters[1] = 0.37e-5;
// fRemovedGammaContributionCorrectionParameters[2] = 0.0003;
- fChargedContributionCorrectionParameters[0] = -0.0231864;
- fChargedContributionCorrectionParameters[1] = 0.112661;
- fChargedContributionCorrectionParameters[2] = 9.2836e-05;
+ fChargedContributionCorrectionParameters[0] = -0.06;
+ fChargedContributionCorrectionParameters[1] = 0.316;
+ fChargedContributionCorrectionParameters[2] = 0.0022;
- fNeutralContributionCorrectionParameters[0] = -9.18740e-03;
- fNeutralContributionCorrectionParameters[1] = 3.58584e-02;
- fNeutralContributionCorrectionParameters[2] = 9.48208e-04;
+ fNeutralContributionCorrectionParameters[0] = -0.003;
+ fNeutralContributionCorrectionParameters[1] = 0.232;
+ fNeutralContributionCorrectionParameters[2] = 0.002;
- fRemovedGammaContributionCorrectionParameters[0] = -2.74323e-03;
- fRemovedGammaContributionCorrectionParameters[1] = 2.71104e-03;
- fRemovedGammaContributionCorrectionParameters[2] = 3.52758e-04;
+ fRemovedGammaContributionCorrectionParameters[0] = 0.001;
+ fRemovedGammaContributionCorrectionParameters[1] = 0.009;
+ fRemovedGammaContributionCorrectionParameters[2] = 0.0;
- fSecondaryContributionCorrectionParameters[0] = 0.075;
- fSecondaryContributionCorrectionParameters[1] = 0.150187;
- fSecondaryContributionCorrectionParameters[2] = 0.00329361;
+ fSecondaryContributionCorrectionParameters[0] = -0.03;
+ fSecondaryContributionCorrectionParameters[1] = 0.221;
+ fSecondaryContributionCorrectionParameters[2] = 0.002;
}
return AliAnalysisEtReconstructed::TrackHitsCalorimeter(track, magField);
}
-Double_t AliAnalysisEtReconstructedPhos::GetChargedContribution(Int_t clusterMult)
-{ // Charged contrib
- if (clusterMult > 0)
- {
- Double_t nPart = fChargedContributionCorrectionParameters[0] + fChargedContributionCorrectionParameters[1]*clusterMult + fChargedContributionCorrectionParameters[2]*clusterMult*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 = 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 = fRemovedGammaContributionCorrectionParameters[0] + fRemovedGammaContributionCorrectionParameters[1]*clusterMult + fRemovedGammaContributionCorrectionParameters[2]*clusterMult*clusterMult;
-
- Double_t contr = nPart*kMEANGAMMA;
-
- return contr;
- }
- return 0;
-}
-
-Double_t AliAnalysisEtReconstructedPhos::GetSecondaryContribution(Int_t clusterMultiplicity)
-{ // Secondary contrib
- if(clusterMultiplicity > 0)
- {
- return fSecondaryContributionCorrectionParameters[0] + fSecondaryContributionCorrectionParameters[1]*clusterMultiplicity + fSecondaryContributionCorrectionParameters[2]*clusterMultiplicity*clusterMultiplicity;
- }
-
- return 0;
-}
-
-
void AliAnalysisEtReconstructedPhos::CreateHistograms()
void CreateHistograms();
- virtual Double_t GetChargedContribution(Int_t clusterMultiplicity);
-
- virtual Double_t GetNeutralContribution(Int_t clusterMultiplicity);
-
- virtual Double_t GetGammaContribution(Int_t clusterMultiplicity);
-
- virtual Double_t GetSecondaryContribution(Int_t clusterMultiplicity);
-
void SetChargedContributionParameters(Double_t par[3])
{
fChargedContributionCorrectionParameters[0] = par[0];
}
return -1;
}
+Bool_t AliAnalysisEtSelector::FromSecondaryInteraction(const TParticle& part, AliStack &stack) const
+{
+ if(0)
+ {
+ // Let's ignore the gamma<->e+/- channel
+ UInt_t pdgCode = TMath::Abs(stack.Particle(part.GetFirstMother())->GetPdgCode());
+
+ if(pdgCode == fgGammaCode || fgEMinusCode)
+ {
+ std::cout << "EM channel" << std::endl;
+ }
+ }
+ Bool_t partVtxSecondary = (
+ TMath::Sqrt(part.Vx()*part.Vx() + part.Vy()*part.Vy()) > fCuts->GetPrimaryVertexCutXY()
+ || TMath::Abs(part.Vz()) > fCuts->GetPrimaryVertexCutZ()
+ )
+ && TMath::Sqrt(part.Vx()*part.Vx()+part.Vy()*part.Vy() + part.Vz()*part.Vz())<(fCuts->GetGeometryPhosDetectorRadius()-10);
+
+ //Let's find suspect decay (typical for secondary interaction)...
+
+ return SuspeciousDecayInChain(211, 111, part, stack);
+
+
+
+
+}
+
+Bool_t AliAnalysisEtSelector::SuspeciousDecayInChain(const Int_t suspectMotherPdg, const Int_t suspectDaughterPdg, const TParticle &part, AliStack& stack) const
+{
+ UInt_t partPdg = TMath::Abs(part.GetPdgCode());
+ if(part.GetFirstMother() == -1)
+ {
+ return kFALSE;
+ }
+ TParticle *mother = stack.Particle(part.GetFirstMother());
+ UInt_t motherPdg = TMath::Abs(mother->GetPdgCode());
+ if((suspectDaughterPdg==partPdg || 2112 == partPdg )&& suspectMotherPdg == motherPdg)
+ {
+ return kTRUE;
+ }
+ return SuspeciousDecayInChain(suspectMotherPdg, suspectDaughterPdg, *mother, stack);
+}
// Cut on geometrical acceptance
virtual Bool_t CutGeometricalAcceptance(const AliVTrack &/*part*/) const { return true; }
+ // From secondary vertex?
+ virtual Bool_t FromSecondaryInteraction(const TParticle& part, AliStack& stack) const;
protected:
TRefArray *fClusterArray; // Array of clusters
AliAnalysisEtCuts *fCuts; // Pointer to the cuts object; DS: also in base class?
-
- Int_t fRunNumber; // run number
+
+ Bool_t SuspeciousDecayInChain(const UInt_t suspectMotherPdg, const UInt_t suspectDaughterPdg, const TParticle& part, AliStack& stack) const;
+
+ Int_t fRunNumber;
Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster) const
{
+
+// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
}
Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const TParticle& part) const
{
+// std::cout << fCuts->GetReconstructedPhosClusterEnergyCut();
return part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
}
TVector3 glVec(gPos);
fGeoUtils->GlobalPos2RelId(glVec, relId);
+ //std::cout << "In phos distance to bad channel cut!" << std::endl;
TVector3 locVec;
fGeoUtils->Global2Local(locVec, glVec, relId[0]);
// std::cout << fGeoUtils << std::endl;
// cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
Int_t nTracksMatched = cluster.GetNTracksMatched();
- if(nTracksMatched == 0) return kFALSE;
+ if(nTracksMatched == 0) return kTRUE;
Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
if(trackMatchedIndex < 0) return kTRUE;
track.Phi() > fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180. &&
track.Phi() < fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
}
+
--- /dev/null
+
+
+#include "AliAnalysisEtTrackMatchCorrections.h"
+
+ClassImp(AliAnalysisEtTrackMatchCorrections);
+
+AliAnalysisEtTrackMatchCorrections::AliAnalysisEtTrackMatchCorrections() : TNamed("TMCorr","TMCorr")
+ ,fChargedContr()
+ ,fNeutralContr()
+ ,fGammaContr()
+ ,fSecondaryContr()
+ ,fMeanCharged(0)
+ ,fMeanNeutral(0)
+ ,fMeanGamma(0)
+ ,fMeanSecondary(0)
+{}
+
+AliAnalysisEtTrackMatchCorrections::AliAnalysisEtTrackMatchCorrections(TString name, TF1 chargedContr, TF1 neutralContr, TF1 gammaContr, TF1 secondaryContr,
+ Double_t meanCharged, Double_t meanNeutral, Double_t meanGammas, Double_t meanSecondary) : TNamed(name,name)
+ ,fChargedContr(chargedContr)
+ ,fNeutralContr(neutralContr)
+ ,fGammaContr(gammaContr)
+ ,fSecondaryContr(secondaryContr)
+ ,fMeanCharged(meanCharged)
+ ,fMeanNeutral(meanNeutral)
+ ,fMeanGamma(meanGammas)
+ ,fMeanSecondary(meanSecondary)
+{}
+
+//! Copy constructor
+AliAnalysisEtTrackMatchCorrections::AliAnalysisEtTrackMatchCorrections(const AliAnalysisEtTrackMatchCorrections &obj) : TNamed(obj)
+ ,fChargedContr(obj.fChargedContr)
+ ,fNeutralContr(obj.fNeutralContr)
+ ,fGammaContr(obj.fGammaContr)
+ ,fSecondaryContr(obj.fSecondaryContr)
+ ,fMeanCharged(obj.fMeanCharged)
+ ,fMeanNeutral(obj.fMeanNeutral)
+ ,fMeanGamma(obj.fMeanGamma)
+ ,fMeanSecondary(obj.fMeanSecondary)
+
+{}
+
+//! Destructor
+AliAnalysisEtTrackMatchCorrections::~AliAnalysisEtTrackMatchCorrections()
+{}
+
+//! Assignment operator
+AliAnalysisEtTrackMatchCorrections& AliAnalysisEtTrackMatchCorrections::operator=(const AliAnalysisEtTrackMatchCorrections &other)
+{
+ if (this != &other)
+ {
+ fChargedContr = other.fChargedContr;
+ fNeutralContr = other.fNeutralContr;
+ fGammaContr = other.fGammaContr;
+ fSecondaryContr = other.fSecondaryContr;
+ fMeanCharged = other.fMeanCharged;
+ fMeanNeutral = other.fMeanNeutral;
+ fMeanGamma = other.fMeanGamma;
+ fMeanSecondary = other.fMeanSecondary;
+
+ }
+ return *this;
+}
+
--- /dev/null
+#ifndef ALIANALYSISETTRACKMATCHCORRECTIONS_H
+#define ALIANALYSISETTRACKMATCHCORRECTIONS_H
+
+
+#include "TNamed.h"
+#include "TF1.h"
+
+class AliAnalysisEtTrackMatchCorrections : public TNamed
+{
+
+public:
+
+//! Default constructor
+ AliAnalysisEtTrackMatchCorrections();
+
+//! Constructor
+ AliAnalysisEtTrackMatchCorrections(TString name, TF1 chargedContr, TF1 neutralContr, TF1 gammaContr, TF1 secondaryContr, Double_t meanCharged, Double_t meanNeutral, Double_t meanGammas, Double_t meanSecondary);
+
+//! Copy constructor
+ AliAnalysisEtTrackMatchCorrections(const AliAnalysisEtTrackMatchCorrections &obj);
+
+//! Destructor
+ virtual ~AliAnalysisEtTrackMatchCorrections();
+
+//! Assignment operator
+ AliAnalysisEtTrackMatchCorrections& operator=(const AliAnalysisEtTrackMatchCorrections& other);
+
+// Getters
+
+ TF1 ChargedContr() const {
+ return fChargedContr;
+ }
+
+ TF1 NeutralContr() const {
+ return fNeutralContr;
+ }
+
+ TF1 GammaContr() const {
+ return fGammaContr;
+ }
+
+ TF1 SecondaryContr() const {
+ return fSecondaryContr;
+ }
+
+ Double_t ChargedContr(int mult) const {
+ return fChargedContr.Eval(mult)*fMeanCharged;
+ }
+
+ Double_t NeutralContr(int mult) const {
+ return fNeutralContr.Eval(mult)*fMeanNeutral;
+ }
+
+ Double_t GammaContr(int mult) const {
+ return -fGammaContr.Eval(mult)*fMeanGamma;
+ }
+
+ Double_t SecondaryContr(int mult) const {
+ return fSecondaryContr.Eval(mult)*fMeanSecondary;
+ }
+
+
+// Setters
+
+ void SetChargedcontr(TF1 chargedContr) {
+ fChargedContr = chargedContr;
+ }
+
+ void SetNeutralcontr(TF1 neutralContr) {
+ fNeutralContr = neutralContr;
+ }
+
+ void SetGammacontr(TF1 gammaContr) {
+ fGammaContr = gammaContr;
+ }
+
+ void SetSecondarycontr(TF1 secondaryContr) {
+ fSecondaryContr = secondaryContr;
+ }
+
+
+private:
+
+ // ChargedContr
+ TF1 fChargedContr;
+ // NeutralContr
+ TF1 fNeutralContr;
+ // GammaContr
+ TF1 fGammaContr;
+ // SecondaryContr
+ TF1 fSecondaryContr;
+
+ // Mean deposited energy from charged particles
+ Double_t fMeanCharged;
+ // Mean deposited energy from neutral particles
+ Double_t fMeanNeutral;
+ // Mean deposited energy from gammas
+ Double_t fMeanGamma;
+ // Mean deposited energy from secondaries
+ Double_t fMeanSecondary;
+
+ // Prohibited
+//! Equality operator
+ bool operator==(const AliAnalysisEtTrackMatchCorrections &other) const;
+ ClassDef(AliAnalysisEtTrackMatchCorrections, 1);
+};
+
+#endif //ALIANALYSISETTRACKMATCHCORRECTIONS_H
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff