#include <TString.h>
#include <TVector2.h>
#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
#include "AliAODVertex.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskEMCALClusterizeFast.h"
#include "AliCDBManager.h"
#include "AliCentrality.h"
-#include "AliEMCALGeoUtils.h"
+#include "AliEMCALDigit.h"
#include "AliEMCALGeometry.h"
+#include "AliEMCALRecPoint.h"
#include "AliEMCALRecoUtils.h"
+#include "AliESDCaloTrigger.h"
#include "AliESDEvent.h"
+#include "AliESDUtils.h"
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
+#include "AliEventplane.h"
#include "AliGeomManager.h"
#include "AliInputEventHandler.h"
#include "AliLog.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
#include "AliMagF.h"
+#include "AliMultiplicity.h"
+#include "AliStack.h"
#include "AliTrackerBase.h"
ClassImp(AliAnalysisTaskEMCALPi0PbPb)
+//________________________________________________________________________
+AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb()
+ : AliAnalysisTaskSE(),
+ fCentVar("V0M"),
+ fCentFrom(0),
+ fCentTo(100),
+ fVtxZMin(-10),
+ fVtxZMax(+10),
+ fUseQualFlag(1),
+ fClusName(),
+ fDoNtuple(0),
+ fDoAfterburner(0),
+ fAsymMax(1),
+ fNminCells(1),
+ fMinE(0.100),
+ fMinErat(0),
+ fMinEcc(0),
+ fGeoName("EMCAL_FIRSTYEARV1"),
+ fMinNClusPerTr(50),
+ fIsoDist(0.2),
+ fTrClassNames(""),
+ fTrCuts(0),
+ fPrimTrCuts(0),
+ fDoTrMatGeom(0),
+ fTrainMode(0),
+ fMarkCells(),
+ fMinL0Time(-1),
+ fMaxL0Time(1024),
+ fMcMode(0),
+ fEmbedMode(0),
+ fGeom(0),
+ fReco(0),
+ fDoPSel(kTRUE),
+ fIsGeoMatsSet(0),
+ fNEvs(0),
+ fOutput(0),
+ fTrClassNamesArr(0),
+ fEsdEv(0),
+ fAodEv(0),
+ fRecPoints(0),
+ fDigits(0),
+ fEsdClusters(0),
+ fEsdCells(0),
+ fAodClusters(0),
+ fAodCells(0),
+ fPtRanges(0),
+ fSelTracks(0),
+ fSelPrimTracks(0),
+ fNAmpInTrigger(0),
+ fAmpInTrigger(0),
+ fNtuple(0),
+ fHeader(0),
+ fPrimVert(0),
+ fSpdVert(0),
+ fTpcVert(0),
+ fClusters(0),
+ fTriggers(0),
+ fMcParts(0),
+ fHCuts(0x0),
+ fHVertexZ(0x0),
+ fHVertexZ2(0x0),
+ fHCent(0x0),
+ fHCentQual(0x0),
+ fHTclsBeforeCuts(0x0),
+ fHTclsAfterCuts(0x0),
+ fHColuRow(0x0),
+ fHColuRowE(0x0),
+ fHCellMult(0x0),
+ fHCellE(0x0),
+ fHCellH(0x0),
+ fHCellM(0x0),
+ fHCellM2(0x0),
+ fHCellFreqNoCut(0x0),
+ fHCellFreqCut100M(0x0),
+ fHCellFreqCut300M(0x0),
+ fHCellFreqE(0x0),
+ fHCellCheckE(0x0),
+ fHClustEccentricity(0),
+ fHClustEtaPhi(0x0),
+ fHClustEnergyPt(0x0),
+ fHClustEnergySigma(0x0),
+ fHClustSigmaSigma(0x0),
+ fHClustNCellEnergyRatio(0x0),
+ fHClustEnergyNCell(0x0),
+ fHPrimTrackPt(0x0),
+ fHPrimTrackEta(0x0),
+ fHPrimTrackPhi(0x0),
+ fHMatchDr(0x0),
+ fHMatchDz(0x0),
+ fHMatchEp(0x0),
+ fHPionEtaPhi(0x0),
+ fHPionMggPt(0x0),
+ fHPionMggAsym(0x0),
+ fHPionMggDgg(0x0)
+{
+ // Constructor.
+}
+
//________________________________________________________________________
AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
: AliAnalysisTaskSE(name),
fMinErat(0),
fMinEcc(0),
fGeoName("EMCAL_FIRSTYEARV1"),
- fMinNClustPerTrack(50),
- fMinPtPerTrack(1.0),
+ fMinNClusPerTr(50),
fIsoDist(0.2),
fTrClassNames(""),
fTrCuts(0),
- fDoTrackMatWithGeom(0),
- fDoConstrain(0),
- fNEvs(0),
+ fPrimTrCuts(0),
+ fDoTrMatGeom(0),
+ fTrainMode(0),
+ fMarkCells(),
+ fMinL0Time(-1),
+ fMaxL0Time(1024),
+ fMcMode(0),
+ fEmbedMode(0),
fGeom(0),
fReco(0),
+ fDoPSel(kTRUE),
+ fIsGeoMatsSet(0),
+ fNEvs(0),
fOutput(0),
fTrClassNamesArr(0),
fEsdEv(0),
fAodEv(0),
fRecPoints(0),
+ fDigits(0),
fEsdClusters(0),
fEsdCells(0),
fAodClusters(0),
fAodCells(0),
fPtRanges(0),
fSelTracks(0),
+ fSelPrimTracks(0),
+ fNAmpInTrigger(0),
+ fAmpInTrigger(0),
fNtuple(0),
fHeader(0),
fPrimVert(0),
fSpdVert(0),
fTpcVert(0),
fClusters(0),
+ fTriggers(0),
+ fMcParts(0),
fHCuts(0x0),
fHVertexZ(0x0),
fHVertexZ2(0x0),
fHClustEnergySigma(0x0),
fHClustSigmaSigma(0x0),
fHClustNCellEnergyRatio(0x0),
+ fHClustEnergyNCell(0x0),
+ fHPrimTrackPt(0x0),
+ fHPrimTrackEta(0x0),
+ fHPrimTrackPhi(0x0),
+ fHMatchDr(0x0),
+ fHMatchDz(0x0),
+ fHMatchEp(0x0),
fHPionEtaPhi(0x0),
fHPionMggPt(0x0),
fHPionMggAsym(0x0),
{
// Constructor.
- if (!name)
- return;
- SetName(name);
- DefineInput(0, TChain::Class());
DefineOutput(1, TList::Class());
- fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks "
+ fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks EMCALTrigger."
"AOD:header,vertices,emcalCells,tracks";
}
delete fReco; fReco = 0;
delete fTrClassNamesArr;
delete fSelTracks;
+ delete fSelPrimTracks;
+ delete [] fAmpInTrigger;
delete [] fHColuRow;
delete [] fHColuRowE;
delete [] fHCellMult;
{
// Create user objects here.
- fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
- fReco = new AliEMCALRecoUtils();
+ cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
+ cout << " fCentVar: " << fCentVar << endl;
+ cout << " fCentFrom: " << fCentFrom << endl;
+ cout << " fCentTo: " << fCentTo << endl;
+ cout << " fVtxZMin: " << fVtxZMin << endl;
+ cout << " fVtxZMax: " << fVtxZMax << endl;
+ cout << " fUseQualFlag: " << fUseQualFlag << endl;
+ cout << " fClusName: \"" << fClusName << "\"" << endl;
+ cout << " fDoNtuple: " << fDoNtuple << endl;
+ cout << " fDoAfterburner: " << fDoAfterburner << endl;
+ cout << " fAsymMax: " << fAsymMax << endl;
+ cout << " fNminCells: " << fNminCells << endl;
+ cout << " fMinE: " << fMinE << endl;
+ cout << " fMinErat: " << fMinErat << endl;
+ cout << " fMinEcc: " << fMinEcc << endl;
+ cout << " fGeoName: \"" << fGeoName << "\"" << endl;
+ cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
+ cout << " fIsoDist: " << fIsoDist << endl;
+ cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
+ cout << " fTrCuts: " << fTrCuts << endl;
+ cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
+ cout << " fDoTrMatGeom: " << fDoTrMatGeom << endl;
+ cout << " fTrainMode: " << fTrainMode << endl;
+ cout << " fMarkCells: " << fMarkCells << endl;
+ cout << " fMinL0Time: " << fMinL0Time << endl;
+ cout << " fMaxL0Time: " << fMaxL0Time << endl;
+ cout << " fMcMode: " << fMcMode << endl;
+ cout << " fEmbedMode: " << fEmbedMode << endl;
+ cout << " fGeom: " << fGeom << endl;
+ cout << " fReco: " << fReco << endl;
+ cout << " fDoPSel: " << fDoPSel << endl;
+
+ if (!fGeom)
+ fGeom = AliEMCALGeometry::GetInstance(fGeoName);
+ else {
+ if (fGeom->GetMatrixForSuperModule(0))
+ fIsGeoMatsSet = kTRUE;
+ }
+ if (!fReco)
+ fReco = new AliEMCALRecoUtils();
fTrClassNamesArr = fTrClassNames.Tokenize(" ");
fOutput = new TList();
- fOutput->SetOwner();
+ fOutput->SetOwner(1);
fSelTracks = new TObjArray;
+ fSelPrimTracks = new TObjArray;
if (fDoNtuple) {
- TFile *f = OpenFile(1);
+ TFile *f = OpenFile(1);
+ TDirectory::TContext context(f);
if (f) {
f->SetCompressionLevel(2);
fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
fNtuple->SetDirectory(f);
- fNtuple->SetAutoFlush(-1024*1024*1024);
- fNtuple->SetAutoSave(-1024*1024*1024);
+ if (fTrainMode) {
+ fNtuple->SetAutoFlush(-2*1024*1024);
+ fNtuple->SetAutoSave(0);
+ } else {
+ fNtuple->SetAutoFlush(-32*1024*1024);
+ fNtuple->SetAutoSave(0);
+ }
+
fHeader = new AliStaHeader;
fNtuple->Branch("header", &fHeader, 16*1024, 99);
fPrimVert = new AliStaVertex;
fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
if (TClass::GetClass("AliStaCluster"))
TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
- fClusters = new TClonesArray("AliStaCluster",1000);
+ fClusters = new TClonesArray("AliStaCluster");
fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
+ if (TClass::GetClass("AliStaTrigger"))
+ TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
+ fTriggers = new TClonesArray("AliStaTrigger");
+ fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
+ if (fMcMode||fEmbedMode) {
+ if (TClass::GetClass("AliStaPart"))
+ TClass::GetClass("AliStaPart")->IgnoreTObjectStreamer();
+ fMcParts = new TClonesArray("AliStaPart");
+ fNtuple->Branch("mcparts", &fMcParts, 8*16*1024, 99);
+ }
}
}
AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
- Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
- Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
+ Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
+ Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
// histograms
TH1::SetDefaultSumw2(kTRUE);
fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
fHVertexZ2->SetXTitle("z [cm]");
fOutput->Add(fHVertexZ2);
- fHCent = new TH1F("hCentBeforeCut","",101,-1,100);
+ fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
fHCent->SetXTitle(fCentVar.Data());
fOutput->Add(fHCent);
- fHCentQual = new TH1F("hCentAfterCut","",101,-1,100);
+ fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
fHCentQual->SetXTitle(fCentVar.Data());
fOutput->Add(fHCentQual);
fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
fOutput->Add(fHCellFreqCut300M[i]);
fOutput->Add(fHCellFreqE[i]);
}
- if (1) {
+ if (!fMarkCells.IsNull()) {
fHCellCheckE = new TH1*[24*48*nsm];
memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
- Int_t tcs[1] = {4102};
- for (UInt_t i = 0; i<sizeof(tcs)/sizeof(Int_t); ++i){
+ TObjArray *cells = fMarkCells.Tokenize(" ");
+ Int_t n = cells->GetEntries();
+ Int_t *tcs = new Int_t[n];
+ for (Int_t i=0;i<n;++i) {
+ TString name(cells->At(i)->GetName());
+ tcs[i]=name.Atoi();
+ }
+ for (Int_t i = 0; i<n; ++i) {
Int_t c=tcs[i];
if (c<24*48*nsm) {
- fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 500, 0, 8);
+ fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
fOutput->Add(fHCellCheckE[i]);
}
}
+ delete cells;
+ delete [] tcs;
}
// histograms for clusters
- fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
- fHClustEccentricity->SetXTitle("#epsilon_{C}");
- fOutput->Add(fHClustEccentricity);
- fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
- fHClustEtaPhi->SetXTitle("#eta");
- fHClustEtaPhi->SetYTitle("#varphi");
- fOutput->Add(fHClustEtaPhi);
- fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
- fHClustEnergyPt->SetXTitle("E [GeV]");
- fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
- fOutput->Add(fHClustEnergyPt);
- fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
- fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
- fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
- fOutput->Add(fHClustEnergySigma);
- fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
- fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
- fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
- fOutput->Add(fHClustSigmaSigma);
- fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
- fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
- fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
- fOutput->Add(fHClustNCellEnergyRatio);
+ if (!fTrainMode) {
+ fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
+ fHClustEccentricity->SetXTitle("#epsilon_{C}");
+ fOutput->Add(fHClustEccentricity);
+ fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
+ fHClustEtaPhi->SetXTitle("#eta");
+ fHClustEtaPhi->SetYTitle("#varphi");
+ fOutput->Add(fHClustEtaPhi);
+ fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
+ fHClustEnergyPt->SetXTitle("E [GeV]");
+ fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
+ fOutput->Add(fHClustEnergyPt);
+ fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
+ fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
+ fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
+ fOutput->Add(fHClustEnergySigma);
+ fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
+ fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
+ fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
+ fOutput->Add(fHClustSigmaSigma);
+ fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
+ fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
+ fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
+ fOutput->Add(fHClustNCellEnergyRatio);
+ fHClustEnergyNCell = new TH2F("hClustEnergyNCell","",200,0,100,50,0,50);
+ fHClustEnergyNCell->SetXTitle("E_{clus}");
+ fHClustEnergyNCell->SetYTitle("N_{cells}");
+ fOutput->Add(fHClustEnergyNCell);
+ }
+
+ // histograms for primary tracks
+ fHPrimTrackPt = new TH1F("hPrimTrackPt",";p_{T} [GeV/c]",500,0,50);
+ fOutput->Add(fHPrimTrackPt);
+ fHPrimTrackEta = new TH1F("hPrimTrackEta",";#eta",40,-2,2);
+ fOutput->Add(fHPrimTrackEta);
+ fHPrimTrackPhi = new TH1F("hPrimTrackPhi",";#varPhi [rad]",63,0,6.3);
+ fOutput->Add(fHPrimTrackPhi);
+ // histograms for track matching
+ fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
+ fOutput->Add(fHMatchDr);
+ fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
+ fOutput->Add(fHMatchDz);
+ fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
+ fOutput->Add(fHMatchEp);
// histograms for pion candidates
- fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
- fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
- fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
- fOutput->Add(fHPionEtaPhi);
- fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
- fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
- fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
- fOutput->Add(fHPionMggPt);
- fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
- fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
- fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
- fOutput->Add(fHPionMggAsym);
- fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
- fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
- fHPionMggDgg->SetYTitle("opening angle [grad]");
- fOutput->Add(fHPionMggDgg);
- const Int_t nbins = 20;
- Double_t xbins[nbins] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12.5,15,20,25,50};
- fPtRanges = new TAxis(nbins-1,xbins);
- for (Int_t i = 0; i<=nbins; ++i) {
- fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
- fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
- if (i==0)
- fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
- else if (i==nbins)
- fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
- else
- fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
- fOutput->Add(fHPionInvMasses[i]);
+ if (!fTrainMode) {
+ fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
+ fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
+ fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
+ fOutput->Add(fHPionEtaPhi);
+ fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
+ fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
+ fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
+ fOutput->Add(fHPionMggPt);
+ fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
+ fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
+ fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
+ fOutput->Add(fHPionMggAsym);
+ fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
+ fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
+ fHPionMggDgg->SetYTitle("opening angle [grad]");
+ fOutput->Add(fHPionMggDgg);
+ const Int_t nbins = 20;
+ Double_t xbins[nbins] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12.5,15,20,25,50};
+ fPtRanges = new TAxis(nbins-1,xbins);
+ for (Int_t i = 0; i<=nbins; ++i) {
+ fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
+ fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
+ if (i==0)
+ fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
+ else if (i==nbins)
+ fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
+ else
+ fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
+ fOutput->Add(fHPionInvMasses[i]);
+ }
}
PostData(1, fOutput);
AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
+ UInt_t offtrigger = 0;
if (fEsdEv) {
am->LoadBranch("AliESDRun.");
am->LoadBranch("AliESDHeader.");
+ UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
+ UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
+ Bool_t desc1 = (mask1 >> 18) & 0x1;
+ Bool_t desc2 = (mask2 >> 18) & 0x1;
+ if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(18)=="EMCAL"
+ AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
+ mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
+ mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
+ return;
+ }
+ offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
} else {
fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
if (!fAodEv) {
return;
}
am->LoadBranch("header");
+ offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
+ }
+ if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
+ AliWarning(Form("EMCAL not in fast only partition"));
+ return;
+ }
+ if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
+ AliWarning("Accessing geometry from OCDB, this is not very efficient!");
+ AliCDBManager *cdb = AliCDBManager::Instance();
+ if (!cdb->IsDefaultStorageSet())
+ cdb->SetDefaultStorage("raw://");
+ Int_t runno = InputEvent()->GetRunNumber();
+ if (runno != cdb->GetRun())
+ cdb->SetRun(runno);
+ AliGeomManager::LoadGeometry();
+ }
+
+ if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
+ Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
+ for (Int_t i=0; i<nsm; ++i) {
+ const TGeoHMatrix *geom = 0;
+ if (fEsdEv)
+ geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
+ else
+ geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
+ if (!geom)
+ continue;
+ geom->Print();
+ fGeom->SetMisalMatrix(geom,i);
+ }
+ fIsGeoMatsSet = kTRUE;
+ }
+
+ if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
+ if (fEsdEv) {
+ const AliESDRun *erun = fEsdEv->GetESDRun();
+ AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
+ erun->GetCurrentDip(),
+ AliMagF::kConvLHC,
+ kFALSE,
+ erun->GetBeamEnergy(),
+ erun->GetBeamType());
+ TGeoGlobalMagField::Instance()->SetField(field);
+ } else {
+ Double_t pol = -1; //polarity
+ Double_t be = -1; //beam energy
+ AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
+ Int_t runno = fAodEv->GetRunNumber();
+ if (runno>=136851 && runno<138275) {
+ pol = -1;
+ be = 2760;
+ btype = AliMagF::kBeamTypeAA;
+ } else if (runno>=138275 && runno<=139517) {
+ pol = +1;
+ be = 2760;
+ btype = AliMagF::kBeamTypeAA;
+ } else {
+ AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
+ }
+ TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
+ }
}
Int_t cut = 1;
fHTclsBeforeCuts->Fill(1+i);
}
- UInt_t res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
- if (res==0)
+ if (fDoPSel && offtrigger==0)
return;
- if (fHCuts->GetBinContent(2)==0) {
- if (fDoTrackMatWithGeom && !AliGeomManager::GetGeometry()) { // get geometry
- AliWarning("Accessing geometry from OCDB, this is not very efficient!");
- AliCDBManager *cdb = AliCDBManager::Instance();
- if (!cdb->IsDefaultStorageSet())
- cdb->SetDefaultStorage("raw://");
- Int_t runno = InputEvent()->GetRunNumber();
- if (runno != cdb->GetRun())
- cdb->SetRun(runno);
- AliGeomManager::LoadGeometry();
- }
-
- if (!fGeom->GetMatrixForSuperModule(0)) { // set misalignment matrices (stored in first event)
- if (fEsdEv) {
- for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i)
- fGeom->SetMisalMatrix(fEsdEv->GetESDRun()->GetEMCALMatrix(i),i);
- } else {
- for (Int_t i=0; i<fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++i)
- fGeom->SetMisalMatrix(fAodEv->GetHeader()->GetEMCALMatrix(i),i);
- }
- }
-
- if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
- if (fEsdEv) {
- const AliESDRun *erun = fEsdEv->GetESDRun();
- AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
- erun->GetCurrentDip(),
- AliMagF::kConvLHC,
- kFALSE,
- erun->GetBeamEnergy(),
- erun->GetBeamType());
- TGeoGlobalMagField::Instance()->SetField(field);
- } else {
- Double_t pol = -1; //polarity
- Double_t be = -1; //beam energy
- AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
- Int_t runno = fAodEv->GetRunNumber();
- if (runno>=136851 && runno<138275) {
- pol = -1;
- be = 2760;
- btype = AliMagF::kBeamTypeAA;
- } else if (runno>=138275 && runno<=139517) {
- pol = +1;
- be = 2760;
- btype = AliMagF::kBeamTypeAA;
- } else {
- AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
- }
- TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
- }
- }
- }
-
fHCuts->Fill(cut++);
const AliCentrality *centP = InputEvent()->GetCentrality();
}
fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
- fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set of if clusters are attached
+ fDigits = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
+ fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
fEsdCells = 0; // will be set if ESD input used
- fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set of if clusters are attached
+ fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
// or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
fAodCells = 0; // will be set if AOD input used
TObjArray *ts = am->GetTasks();
cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
if (cltask && cltask->GetClusters()) {
- fRecPoints = const_cast<TObjArray*>(cltask->GetClusters());
+ fRecPoints = cltask->GetClusters();
+ fDigits = cltask->GetDigits();
clusattached = cltask->GetAttachClusters();
if (cltask->GetCalibData()!=0)
recalibrated = kTRUE;
}
}
- if (1 && AODEvent() && !fClusName.IsNull()) {
- TList *l = AODEvent()->GetList();
- TClonesArray *clus = 0;
+ if (1 && !fClusName.IsNull()) {
+ TList *l = 0;
+ if (AODEvent())
+ l = AODEvent()->GetList();
+ else if (fAodEv)
+ l = fAodEv->GetList();
if (l) {
- clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
- fAodClusters = clus;
+ fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
}
}
if (!clusattached)
am->LoadBranch("CaloClusters");
TList *l = fEsdEv->GetList();
- TClonesArray *clus = 0;
if (l) {
- clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
- fEsdClusters = clus;
+ fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
}
}
if (1) {
if (!clusattached)
am->LoadBranch("caloClusters");
TList *l = fAodEv->GetList();
- TClonesArray *clus = 0;
if (l) {
- clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
- fAodClusters = clus;
+ fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
}
}
if (1) {
if (1) {
AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
+ AliDebug(2,Form("fDigits set: %p", fDigits));
AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
if (fDoAfterburner)
ClusterAfterburner();
+ CalcMcInfo();
+ CalcCaloTriggers();
+ CalcPrimTracks();
CalcTracks();
CalcClusterProps();
FillCellHists();
- FillClusHists();
- FillPionHists();
- FillOtherHists();
+ if (!fTrainMode) {
+ FillClusHists();
+ FillPionHists();
+ FillOtherHists();
+ }
+ FillMcHists();
FillNtuple();
+ if (fTrainMode) {
+ fSelTracks->Clear();
+ fSelPrimTracks->Clear();
+ if (fMcParts)
+ fMcParts->Clear();
+ if (fTriggers)
+ fTriggers->Clear();
+ if (fClusters)
+ fClusters->Clear();
+ }
+
PostData(1, fOutput);
}
if (fNtuple) {
TFile *f = OpenFile(1);
+ TDirectory::TContext context(f);
if (f)
fNtuple->Write();
}
- AliInfo(Form("%s: Accepted %lld events", GetName(), fNEvs));
+ AliInfo(Form("%s: Accepted %lld events ", GetName(), fNEvs));
}
//________________________________________________________________________
-void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
+void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
+{
+ // Calculate triggers
+
+ if (fAodEv)
+ return; // information not available in AOD
+
+ if (!fTriggers)
+ return;
+
+ fTriggers->Clear();
+
+ AliVCaloCells *cells = fEsdCells;
+ if (!cells)
+ cells = fAodCells;
+ if (!cells)
+ return;
+
+ Int_t ncells = cells->GetNumberOfCells();
+ if (ncells<=0)
+ return;
+
+ if (ncells>fNAmpInTrigger) {
+ delete [] fAmpInTrigger;
+ fAmpInTrigger = new Float_t[ncells];
+ fNAmpInTrigger = ncells;
+ }
+ for (Int_t i=0;i<ncells;++i)
+ fAmpInTrigger[i] = 0;
+
+ std::map<Short_t,Short_t> map;
+ for (Short_t pos=0;pos<ncells;++pos) {
+ Short_t id = cells->GetCellNumber(pos);
+ map[id]=pos;
+ }
+
+ AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+ am->LoadBranch("EMCALTrigger.");
+
+ AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
+ if (!triggers)
+ return;
+ if (triggers->GetEntries()<=0)
+ return;
+
+ triggers->Reset();
+ Int_t ntrigs=0;
+ while (triggers->Next()) {
+ Int_t gCol=0, gRow=0, ntimes=0;
+ triggers->GetPosition(gCol,gRow);
+ triggers->GetNL0Times(ntimes);
+ if (ntimes<1)
+ continue;
+ Float_t amp=0;
+ triggers->GetAmplitude(amp);
+ Int_t find = -1;
+ fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
+ if (find<0)
+ continue;
+ Int_t cidx[4] = {-1};
+ Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
+ if (!ret)
+ continue;
+ Int_t trgtimes[25];
+ triggers->GetL0Times(trgtimes);
+ Int_t mintime = trgtimes[0];
+ Int_t maxtime = trgtimes[0];
+ Bool_t trigInTimeWindow = 0;
+ for (Int_t i=0;i<ntimes;++i) {
+ if (trgtimes[i]<mintime)
+ mintime = trgtimes[i];
+ if (maxtime<trgtimes[i])
+ maxtime = trgtimes[i];
+ if ((fMinL0Time<=trgtimes[i]) && (fMaxL0Time>=trgtimes[i]))
+ trigInTimeWindow = 1;
+ }
+
+ Double_t tenergy = 0;
+ Double_t tphi=0;
+ Double_t teta=0;
+ for (Int_t i=0;i<3;++i) {
+ Short_t pos = -1;
+ std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
+ if (it!=map.end())
+ pos = it->second;
+ if (pos<0)
+ continue;
+ if (trigInTimeWindow)
+ fAmpInTrigger[pos] = amp;
+ Float_t eta=-1, phi=-1;
+ fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
+ Double_t en= cells->GetAmplitude(pos);
+ tenergy+=en;
+ teta+=eta*en;
+ tphi+=phi*en;
+ }
+
+ if (tenergy<=0)
+ continue;
+
+ teta/=tenergy;
+ tphi/=tenergy;
+
+ AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
+ trignew->fE = tenergy;
+ trignew->fEta = teta;
+ trignew->fPhi = tphi;
+ trignew->fAmp = amp;
+ trignew->fMinTime = mintime;
+ trignew->fMaxTime = maxtime;
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
+{
+ // Calculate cluster properties
+
+ if (!fClusters)
+ return;
+
+ fClusters->Clear();
+
+ AliVCaloCells *cells = fEsdCells;
+ if (!cells)
+ cells = fAodCells;
+ if (!cells)
+ return;
+
+ TObjArray *clusters = fEsdClusters;
+ if (!clusters)
+ clusters = fAodClusters;
+ if (!clusters)
+ return;
+
+ Int_t ncells = cells->GetNumberOfCells();
+ Int_t nclus = clusters->GetEntries();
+ Int_t ntrks = fSelTracks->GetEntries();
+ Bool_t btracks[6][ntrks];
+ memset(btracks,0,sizeof(btracks));//todo
+
+ std::map<Short_t,Short_t> map;
+ for (Short_t pos=0;pos<ncells;++pos) {
+ Short_t id = cells->GetCellNumber(pos);
+ map[id]=pos;
+ }
+
+ TObjArray filtMcParts;
+ if (fMcParts) {
+ Int_t nmc = fMcParts->GetEntries();
+ for (Int_t i=0; i<nmc; ++i) {
+ AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
+ if (pa->OnEmcal())
+ filtMcParts.Add(pa);
+ }
+ }
+
+ for(Int_t i=0, ncl=0; i<nclus; ++i) {
+ AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
+
+ if (!clus)
+ continue;
+ if (!clus->IsEMCAL())
+ continue;
+ if (clus->E()<fMinE)
+ continue;
+
+ Float_t clsPos[3] = {0};
+ clus->GetPosition(clsPos);
+ TVector3 clsVec(clsPos);
+ Double_t vertex[3] = {0};
+ InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+ TLorentzVector clusterVec;
+ clus->GetMomentum(clusterVec,vertex);
+ Double_t clsEta = clusterVec.Eta();
+
+ AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
+ cl->fE = clus->E();
+ cl->fR = clsVec.Perp();
+ cl->fEta = clsVec.Eta();
+ cl->fPhi = clsVec.Phi();
+ cl->fN = clus->GetNCells();
+ cl->fN1 = GetNCells(clus,0.100);
+ cl->fN3 = GetNCells(clus,0.300);
+ Short_t id = -1;
+ Double_t emax = GetMaxCellEnergy(clus, id);
+ cl->fIdMax = id;
+ cl->fEmax = emax;
+ cl->fTmax = cells->GetCellTime(id);
+ if (clus->GetDistanceToBadChannel()<10000)
+ cl->fDbc = clus->GetDistanceToBadChannel();
+ if (!TMath::IsNaN(clus->GetDispersion()))
+ cl->fDisp = clus->GetDispersion();
+ if (!TMath::IsNaN(clus->GetM20()))
+ cl->fM20 = clus->GetM20();
+ if (!TMath::IsNaN(clus->GetM02()))
+ cl->fM02 = clus->GetM02();
+ Double_t maxAxis = 0, minAxis = 0;
+ GetSigma(clus,maxAxis,minAxis);
+ clus->SetTOF(maxAxis); // store sigma in TOF
+ cl->fSig = maxAxis;
+ Double_t clusterEcc = 0;
+ if (maxAxis > 0)
+ clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
+ clus->SetChi2(clusterEcc); // store ecc in chi2
+ cl->fEcc = clusterEcc;
+ cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
+ cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
+ cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
+ cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
+ cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
+
+ if (fAmpInTrigger) { // fill trigger info if present
+ Double_t trigpen = 0;
+ Double_t trignen = 0;
+ for(Int_t j=0; j<cl->fN; ++j) {
+ Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
+ Short_t pos = -1;
+ std::map<Short_t,Short_t>::iterator it = map.find(cid);
+ if (it!=map.end())
+ pos = it->second;
+ if (pos<0)
+ continue;
+ if (fAmpInTrigger[pos]>0)
+ trigpen += cells->GetAmplitude(pos);
+ else if (fAmpInTrigger[pos]<0)
+ trignen += cells->GetAmplitude(pos);
+ }
+ if (trigpen>0) {
+ cl->fIsTrigM = 1;
+ cl->fTrigE = trigpen;
+ }
+ if (trignen>0) {
+ cl->fIsTrigM = 1;
+ cl->fTrigMaskE = trignen;
+ }
+ }
+ cl->fIsShared = IsShared(clus);
+
+ // track matching
+ Double_t mind2 = 1e10;
+ for(Int_t j = 0; j<ntrks; ++j) {
+ AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
+ if (!track)
+ continue;
+
+ if (TMath::Abs(clsEta-track->Eta())>0.5)
+ continue;
+
+ TVector3 vec(clsPos);
+ Int_t index = (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
+ if (btracks[index-4][j]) {
+ continue;
+ }
+
+ Float_t tmpR=-1, tmpZ=-1;
+ if (!fDoTrMatGeom) {
+ AliExternalTrackParam *tParam = 0;
+ if (fEsdEv) {
+ AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
+ tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
+ } else
+ tParam = new AliExternalTrackParam(track);
+
+ Double_t bfield[3] = {0};
+ track->GetBxByBz(bfield);
+ Double_t alpha = (index+0.5)*20*TMath::DegToRad();
+ vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
+ tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
+ Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
+ if (!ret) {
+ btracks[index-4][j]=1;
+ delete tParam;
+ continue;
+ }
+ Double_t trkPos[3] = {0};
+ tParam->GetXYZ(trkPos); //Get the extrapolated global position
+ tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2) +
+ TMath::Power(clsPos[1]-trkPos[1],2) +
+ TMath::Power(clsPos[2]-trkPos[2],2) );
+ tmpZ = clsPos[2]-trkPos[2];
+ delete tParam;
+ } else {
+ if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
+ continue;
+ AliExternalTrackParam tParam(track);
+ if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
+ continue;
+ }
+
+ Double_t d2 = tmpR;
+ if (mind2>d2) {
+ mind2=d2;
+ cl->fTrDz = tmpZ;
+ cl->fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
+ cl->fTrEp = clus->E()/track->P();
+ cl->fIsTrackM = 1;
+ }
+ }
+
+ if (cl->fIsTrackM) {
+ fHMatchDr->Fill(cl->fTrDr);
+ fHMatchDz->Fill(cl->fTrDz);
+ fHMatchEp->Fill(cl->fTrEp);
+ }
+
+ //mc matching
+ if (fMcParts) {
+ Int_t nmc = filtMcParts.GetEntries();
+ Double_t diffR2 = 1e9;
+ AliStaPart *msta = 0;
+ for (Int_t j=0; j<nmc; ++j) {
+ AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
+ Double_t t1=clsVec.Eta()-pa->fVEta;
+ Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
+ Double_t tmp = t1*t1+t2*t2;
+ if (tmp<diffR2) {
+ diffR2 = tmp;
+ msta = pa;
+ }
+ }
+ if (diffR2<10 && msta!=0) {
+ cl->fMcLabel = msta->fLab;
+ }
+ }
+
+ cl->fEmbE = 0;
+ if (fDigits && fEmbedMode) {
+ for(Int_t j=0; j<cl->fN; ++j) {
+ Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
+ Short_t pos = -1;
+ std::map<Short_t,Short_t>::iterator it = map.find(cid);
+ if (it!=map.end())
+ pos = it->second;
+ if (pos<0)
+ continue;
+ AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
+ if (!digit)
+ continue;
+ if (digit->GetId() != cid) {
+ AliError(Form("Ids should be equal: %d %d", cid, digit->GetId()));
+ continue;
+ }
+ if (digit->GetType()<-1) {
+ cl->fEmbE += digit->GetChi2();
+ }
+ }
+ }
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
{
// Calculate track properties.
- fSelTracks->Clear();
+ fSelPrimTracks->Clear();
AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
TClonesArray *tracks = 0;
return;
AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
- Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-0.25;
- Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+0.25;
+ Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
+ Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
if (fEsdEv) {
- if (fDoConstrain)
- fSelTracks->SetOwner(kTRUE);
+ fSelPrimTracks->SetOwner(kTRUE);
am->LoadBranch("PrimaryVertex.");
const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
am->LoadBranch("SPDVertex.");
continue;
if (track->Phi()<phimin||track->Phi()>phimax)
continue;
- if(track->GetTPCNcls()<fMinNClustPerTrack)
- continue;
-
- if (!fDoConstrain) {
- fSelTracks->Add(track);
- continue;
- }
AliESDtrack copyt(*track);
Double_t bfield[3];
aTrack->SetChi2perNDF(-1);
aTrack->SetFlags(copyt.GetStatus());
aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
- fSelTracks->Add(aTrack);
+ fSelPrimTracks->Add(aTrack);
}
} else {
Int_t ntracks = tracks->GetEntries();
continue;
if (track->Phi()<phimin||track->Phi()>phimax)
continue;
- if(track->GetTPCNcls()<fMinNClustPerTrack)
+ if(track->GetTPCNcls()<fMinNClusPerTr)
continue;
-
- if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
- AliExternalTrackParam tParam(track);
- if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
- Double_t trkPos[3];
- tParam.GetXYZ(trkPos);
- track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
- }
- }
- fSelTracks->Add(track);
+ //todo: Learn how to set/filter AODs for prim/sec tracks
+ fSelPrimTracks->Add(track);
}
}
}
//________________________________________________________________________
-void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
+void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
{
- // Calculate cluster properties
-
- TObjArray *clusters = fEsdClusters;
- if (!clusters)
- clusters = fAodClusters;
- if (!clusters)
- return;
-
- Int_t nclus = clusters->GetEntries();
- Int_t ntrks = fSelTracks->GetEntries();
+ // Calculate track properties (including secondaries).
- for(Int_t i = 0; i<nclus; ++i) {
- fClusProps[i].Reset();
+ fSelTracks->Clear();
- AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
- if (!clus)
- continue;
- if (!clus->IsEMCAL())
- continue;
- if (clus->E()<fMinE)
- continue;
+ AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+ TClonesArray *tracks = 0;
+ if (fEsdEv) {
+ am->LoadBranch("Tracks");
+ TList *l = fEsdEv->GetList();
+ tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
+ } else {
+ am->LoadBranch("tracks");
+ TList *l = fAodEv->GetList();
+ tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
+ }
- Float_t clsPos[3] = {0};
- clus->GetPosition(clsPos);
- TVector3 clsVec(clsPos);
- Double_t vertex[3] = {0};
- InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
- TLorentzVector clusterVec;
- clus->GetMomentum(clusterVec,vertex);
- Double_t clsEta = clusterVec.Eta();
+ if (!tracks)
+ return;
- Double_t mind2 = 1e10;
- for(Int_t j = 0; j<ntrks; ++j) {
- AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
+ if (fEsdEv) {
+ const Int_t Ntracks = tracks->GetEntries();
+ for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
+ AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
+ if (!track) {
+ AliWarning(Form("Could not receive track %d\n", iTracks));
+ continue;
+ }
+ if (fTrCuts && !fTrCuts->IsSelected(track))
+ continue;
+ Double_t eta = track->Eta();
+ if (eta<-1||eta>1)
+ continue;
+ fSelTracks->Add(track);
+ }
+ } else {
+ Int_t ntracks = tracks->GetEntries();
+ for (Int_t i=0; i<ntracks; ++i) {
+ AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
if (!track)
continue;
- Double_t pt = track->Pt();
- if (pt<fMinPtPerTrack)
+ Double_t eta = track->Eta();
+ if (eta<-1||eta>1)
continue;
- if (TMath::Abs(clsEta-track->Eta())>0.5)
+ if(track->GetTPCNcls()<fMinNClusPerTr)
continue;
- Float_t tmpR=-1, tmpZ=-1;
- if (!fDoTrackMatWithGeom) {
-
- AliExternalTrackParam *tParam = 0;
- if (fEsdEv) {
- AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
- tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
- } else
- tParam = new AliExternalTrackParam(track);
-
- Double_t bfield[3];
- track->GetBxByBz(bfield);
- TVector3 vec(clsPos);
- Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
- vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
- tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
- Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
- if (!ret) {
- delete tParam;
- continue;
- }
- Double_t trkPos[3];
- tParam->GetXYZ(trkPos); //Get the extrapolated global position
- tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
- tmpZ = clsPos[2]-trkPos[2];
- delete tParam;
- } else {
- if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
- continue;
+ if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
AliExternalTrackParam tParam(track);
- if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
- continue;
- }
-
- Double_t d2 = tmpR;
- if (mind2>d2) {
- mind2=d2;
- fClusProps[i].fTrIndex = j;
- fClusProps[i].fTrDz = tmpZ;
- fClusProps[i].fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
- fClusProps[i].fTrDist = d2;
- fClusProps[i].fTrEp = clus->E()/track->P();
+ if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
+ Double_t trkPos[3];
+ tParam.GetXYZ(trkPos);
+ track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
+ }
}
+ fSelTracks->Add(track);
}
-
- if (0 && (fClusProps[i].fTrIndex>=0)) {
- cout << i << " " << fClusProps[i].fTrIndex << ": Dr " << fClusProps[i].fTrDr << " " << " Dz " << fClusProps[i].fTrDz << endl;
- }
-
- fClusProps[i].fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
- fClusProps[i].fTrLowPtIso = 0;
- fClusProps[i].fCellIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
}
}
//________________________________________________________________________
-void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
+void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
{
// Run custer reconstruction afterburner.
continue;
TLorentzVector clusterVec;
clus->GetMomentum(clusterVec,vertex);
- Double_t maxAxis = 0, minAxis = 0;
- GetSigma(clus,maxAxis,minAxis);
- clus->SetTOF(maxAxis); // store sigma in TOF
- Double_t clusterEcc = 0;
- if (maxAxis > 0)
- clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
- clus->SetChi2(clusterEcc); // store ecc in chi2
+ Double_t maxAxis = clus->GetTOF(); //sigma
+ Double_t clusterEcc = clus->Chi2(); //eccentricity
fHClustEccentricity->Fill(clusterEcc);
fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
+ //====
+ fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
+{
+ // Get Mc truth particle information.
+
+ if (!fMcParts)
+ return;
+
+ fMcParts->Clear();
+
+ AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
+ Double_t etamin = -0.7;
+ Double_t etamax = +0.7;
+ Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
+ Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
+
+ if (fAodEv) {
+ AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+ am->LoadBranch(AliAODMCParticle::StdBranchName());
+ TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!tca)
+ return;
+
+ Int_t nents = tca->GetEntries();
+ for(int it=0; it<nents; ++it) {
+ AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
+ part->Print();
+
+ // pion or eta meson or direct photon
+ if(part->GetPdgCode() == 111) {
+ } else if(part->GetPdgCode() == 221) {
+ } else if(part->GetPdgCode() == 22 ) {
+ } else
+ continue;
+
+ // primary particle
+ Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
+ if(dR > 1.0)
+ continue;
+
+ // kinematic cuts
+ Double_t pt = part->Pt() ;
+ if (pt<0.5)
+ continue;
+ Double_t eta = part->Eta();
+ if (eta<etamin||eta>etamax)
+ continue;
+ Double_t phi = part->Phi();
+ if (phi<phimin||phi>phimax)
+ continue;
+
+ ProcessDaughters(part, it, tca);
+ }
+ return;
+ }
+
+ AliMCEvent *mcEvent = MCEvent();
+ if (!mcEvent)
+ return;
+
+ const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
+ if (!evtVtx)
+ return;
+
+ mcEvent->PreReadAll();
+
+ Int_t nTracks = mcEvent->GetNumberOfPrimaries();
+ for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
+ AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
+ if (!mcP)
+ continue;
+
+ // pion or eta meson or direct photon
+ if(mcP->PdgCode() == 111) {
+ } else if(mcP->PdgCode() == 221) {
+ } else if(mcP->PdgCode() == 22 ) {
+ } else
+ continue;
+
+ // primary particle
+ if(mcP->GetMother()>=0 && mcP->GetMother()<nTracks)
+ continue;
+ Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
+ (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
+ if(dR > 1.0)
+ continue;
+
+ // kinematic cuts
+ Double_t pt = mcP->Pt() ;
+ if (pt<0.5)
+ continue;
+ Double_t eta = mcP->Eta();
+ if (eta<etamin||eta>etamax)
+ continue;
+ Double_t phi = mcP->Phi();
+ if (phi<phimin||phi>phimax)
+ continue;
+
+ ProcessDaughters(mcP, iTrack, mcEvent);
}
}
fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
+ Float_t v0CorrR = 0;
+ fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
+ const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
+ if (mult)
+ fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
+ fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
}
AliCentrality *cent = InputEvent()->GetCentrality();
fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
fHeader->fCqual = cent->GetQuality();
-
+
+ AliEventplane *ep = InputEvent()->GetEventplane();
+ if (ep) {
+ if (ep->GetQVector())
+ fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
+ else
+ fHeader->fPsi = -1;
+ if (ep->GetQsub1()&&ep->GetQsub2())
+ fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
+ else
+ fHeader->fPsiRes = 0;
+ }
+
Double_t val = 0;
TString trgclasses(fHeader->fFiredTriggers);
for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
}
fHeader->fTcls = (UInt_t)val;
+ fHeader->fNSelTr = fSelTracks->GetEntries();
+ fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
+ fHeader->fNSelPrimTr1 = 0;
+ fHeader->fNSelPrimTr2 = 0;
+ for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
+ AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
+ if(track->Pt()>1)
+ ++fHeader->fNSelPrimTr1;
+ if(track->Pt()>2)
+ ++fHeader->fNSelPrimTr2;
+ }
+
+ fHeader->fNCells = 0;
+ fHeader->fNCells1 = 0;
+ fHeader->fNCells2 = 0;
+ fHeader->fNCells5 = 0;
+ fHeader->fMaxCellE = 0;
+
+ AliVCaloCells *cells = fEsdCells;
+ if (!cells)
+ cells = fAodCells;
+
+ if (cells) {
+ Int_t ncells = cells->GetNumberOfCells();
+ for(Int_t j=0; j<ncells; ++j) {
+ Double_t cellen = cells->GetAmplitude(j);
+ if (cellen>1)
+ ++fHeader->fNCells1;
+ if (cellen>2)
+ ++fHeader->fNCells2;
+ if (cellen>5)
+ ++fHeader->fNCells5;
+ if (cellen>fHeader->fMaxCellE)
+ fHeader->fMaxCellE = cellen;
+ }
+ fHeader->fNCells = ncells;
+ }
+
+ fHeader->fNClus = 0;
+ fHeader->fNClus1 = 0;
+ fHeader->fNClus2 = 0;
+ fHeader->fNClus5 = 0;
+ fHeader->fMaxClusE = 0;
+
+ TObjArray *clusters = fEsdClusters;
+ if (!clusters)
+ clusters = fAodClusters;
+
+ if (clusters) {
+ Int_t nclus = clusters->GetEntries();
+ for(Int_t j=0; j<nclus; ++j) {
+ AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
+ if (!clus->IsEMCAL())
+ continue;
+ Double_t clusen = clus->E();
+ if (clusen>1)
+ ++fHeader->fNClus1;
+ if (clusen>2)
+ ++fHeader->fNClus2;
+ if (clusen>5)
+ ++fHeader->fNClus5;
+ if (clusen>fHeader->fMaxClusE)
+ fHeader->fMaxClusE = clusen;
+ }
+ fHeader->fNClus = nclus;
+ }
+
if (fAodEv) {
am->LoadBranch("vertices");
AliAODVertex *pv = fAodEv->GetPrimaryVertex();
FillVertex(fTpcVert, tv);
}
- TObjArray *clusters = fEsdClusters;
- if (!clusters)
- clusters = fAodClusters;
- if (!clusters)
- return;
-
- fClusters->Clear();
- Int_t nclus = clusters->GetEntries();
- for(Int_t i=0,ncl=0; i<nclus; ++i) {
- AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
- if (!clus)
- continue;
- if (!clus->IsEMCAL())
- continue;
- if (clus->E()<fMinE)
- continue;
-
- AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
- cl->fE = clus->E();
- Float_t pos[3];
- clus->GetPosition(pos);
- TVector3 vpos(pos);
- cl->fR = vpos.Perp();
- cl->fEta = vpos.Eta();
- cl->fPhi = vpos.Phi();
- cl->fN = clus->GetNCells();
- cl->fN1 = GetNCells(clus,0.100);
- cl->fN3 = GetNCells(clus,0.300);
- Short_t id = -1;
- Double_t emax = GetMaxCellEnergy(clus, id);
- cl->fIdMax = id;
- cl->fEmax = emax;
- cl->fDbc = clus->GetDistanceToBadChannel();;
- cl->fDisp = clus->GetDispersion();
- cl->fM20 = clus->GetM20();
- cl->fM02 = clus->GetM02();
- cl->fEcc = clus->Chi2(); //eccentricity
- cl->fSig = clus->GetTOF(); //sigma
- cl->fTrDz = fClusProps[i].fTrDz;
- cl->fTrDr = fClusProps[i].fTrDr;;
- cl->fTrEp = fClusProps[i].fTrEp;;
- cl->fTrIso = fClusProps[i].fTrIso;
- cl->fCeIso = fClusProps[i].fCellIso;
- }
fNtuple->Fill();
}
}
}
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
+{
+ // Fill additional MC information histograms.
+
+ if (!fMcParts)
+ return;
+
+ // check if aod or esd mc mode and the fill histos
+}
+
//________________________________________________________________________
void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
{
// Fill other histograms.
+
+ for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
+ AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
+ if(!track)
+ continue;
+ fHPrimTrackPt->Fill(track->Pt());
+ fHPrimTrackEta->Fill(track->Eta());
+ fHPrimTrackPhi->Fill(track->Phi());
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
+{
+ // Fill track histograms.
+
+ if (fSelPrimTracks) {
+ for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
+ AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
+ if(!track)
+ continue;
+ fHPrimTrackPt->Fill(track->Pt());
+ fHPrimTrackEta->Fill(track->Eta());
+ fHPrimTrackPhi->Fill(track->Phi());
+ }
+ }
}
//__________________________________________________________________________________________________
Int_t ncells = cells->GetNumberOfCells();
for (Int_t i = 0; i<ncells; ++i) {
Int_t absID = TMath::Abs(cells->GetCellNumber(i));
- Double_t cellE = cells->GetAmplitude(i);
- Float_t eta, phi;
+ Float_t eta=-1, phi=-1;
fGeom->EtaPhiFromIndex(absID,eta,phi);
- Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
+ Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
if(dist>rad2)
continue;
+ Double_t cellE = cells->GetAmplitude(i);
cellIsolation += cellE;
}
return cellIsolation;
}
//________________________________________________________________________
-Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(AliVCluster *cluster, Short_t &id) const
+Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
+{
+ // Get maximum energy of attached cell.
+
+ Double_t ret = 0;
+ Int_t ncells = cluster->GetNCells();
+ if (fEsdCells) {
+ for (Int_t i=0; i<ncells; i++) {
+ Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
+ ret += e;
+ }
+ } else {
+ for (Int_t i=0; i<ncells; i++) {
+ Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
+ ret += e;
+ }
+ }
+ return ret;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
{
// Get maximum energy of attached cell.
}
//________________________________________________________________________
-void AliAnalysisTaskEMCALPi0PbPb::GetSigma(AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
+void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
{
// Calculate the (E) weighted variance along the longer (eigen) axis.
if (ncells==1)
return;
- TVector3 pos;
for(Int_t j=0; j<ncells; ++j) {
Int_t id = TMath::Abs(c->GetCellAbsId(j));
- fGeom->GetGlobal(id,pos);
Double_t cellen = cells->GetCellAmplitude(id);
+ TVector3 pos;
+ fGeom->GetGlobal(id,pos);
Xc += cellen*pos.X();
Yc += cellen*pos.Y();
Sxx += cellen*pos.X()*pos.X();
}
//________________________________________________________________________
-Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(AliVCluster *c, Double_t emin) const
+Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
{
// Calculate number of attached cells above emin.
}
//________________________________________________________________________
-Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
+Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
{
// Compute isolation based on tracks.
Double_t trkIsolation = 0;
Double_t rad2 = radius*radius;
- Int_t ntrks = fSelTracks->GetEntries();
+ Int_t ntrks = fSelPrimTracks->GetEntries();
for(Int_t j = 0; j<ntrks; ++j) {
AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
if (!track)
continue;
+ if (track->Pt()<pt)
+ continue;
Float_t eta = track->Eta();
Float_t phi = track->Phi();
- Double_t phidiff = TVector2::Phi_0_2pi(phi-cPhi);
+ Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
if(dist>rad2)
continue;
}
return trkIsolation;
}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
+{
+ // Returns if cluster shared across super modules.
+
+ AliVCaloCells *cells = fEsdCells;
+ if (!cells)
+ cells = fAodCells;
+ if (!cells)
+ return 0;
+
+ Int_t n = -1;
+ Int_t ncells = c->GetNCells();
+ for(Int_t j=0; j<ncells; ++j) {
+ Int_t id = TMath::Abs(c->GetCellAbsId(j));
+ Int_t got = id / (24*48);
+ if (n==-1) {
+ n = got;
+ continue;
+ }
+ if (got!=n)
+ return 1;
+ }
+ return 0;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
+{
+ // Print recursively daughter information.
+
+ if (!p || !arr)
+ return;
+
+ const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
+ if (!amc)
+ return;
+ for (Int_t i=0; i<level; ++i) printf(" ");
+ amc->Print();
+
+ Int_t n = amc->GetNDaughters();
+ for (Int_t i=0; i<n; ++i) {
+ Int_t d = amc->GetDaughter(i);
+ const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
+ PrintDaughters(dmc,arr,level+1);
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
+{
+ // Print recursively daughter information.
+
+ if (!p || !arr)
+ return;
+
+ for (Int_t i=0; i<level; ++i) printf(" ");
+ Int_t d1 = p->GetFirstDaughter();
+ Int_t d2 = p->GetLastDaughter();
+ printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
+ p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
+ if (d1<0)
+ return;
+ if (d2<0)
+ d2=d1;
+ for (Int_t i=d1;i<=d2;++i) {
+ const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
+ PrintDaughters(dmc,arr,level+1);
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
+{
+ // Print track ref array.
+
+ if (!p)
+ return;
+
+ Int_t n = p->GetNumberOfTrackReferences();
+ for (Int_t i=0; i<n; ++i) {
+ AliTrackReference *ref = p->GetTrackReference(i);
+ if (!ref)
+ continue;
+ ref->SetUserId(ref->DetectorId());
+ ref->Print();
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
+{
+ // Process and create daughters.
+
+ if (!p || !arr)
+ return;
+
+ AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
+ if (!amc)
+ return;
+
+ //amc->Print();
+
+ Int_t nparts = arr->GetEntries();
+ Int_t nents = fMcParts->GetEntries();
+
+ AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
+ newp->fPt = amc->Pt();
+ newp->fEta = amc->Eta();
+ newp->fPhi = amc->Phi();
+ if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
+ TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
+ newp->fVR = vec.Perp();
+ newp->fVEta = vec.Eta();
+ newp->fVPhi = vec.Phi();
+ }
+ newp->fPid = amc->PdgCode();
+ newp->fLab = nents;
+ Int_t moi = amc->GetMother();
+ if (moi>=0&&moi<nparts) {
+ const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
+ moi = mmc->GetUniqueID();
+ }
+ newp->fMo = moi;
+ p->SetUniqueID(nents);
+
+ // TODO: Determine which detector was hit
+ //newp->fDet = ???
+
+ Int_t n = amc->GetNDaughters();
+ for (Int_t i=0; i<n; ++i) {
+ Int_t d = amc->GetDaughter(i);
+ if (d<=index || d>=nparts)
+ continue;
+ AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
+ ProcessDaughters(dmc,d,arr);
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
+{
+ // Process and create daughters.
+
+ if (!p || !arr)
+ return;
+
+ Int_t d1 = p->GetFirstDaughter();
+ Int_t d2 = p->GetLastDaughter();
+ if (0) {
+ printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
+ index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
+ PrintTrackRefs(p);
+ }
+ Int_t nents = fMcParts->GetEntries();
+
+ AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
+ newp->fPt = p->Pt();
+ newp->fEta = p->Eta();
+ newp->fPhi = p->Phi();
+ if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
+ TVector3 vec(p->Xv(),p->Yv(),p->Zv());
+ newp->fVR = vec.Perp();
+ newp->fVEta = vec.Eta();
+ newp->fVPhi = vec.Phi();
+ }
+ newp->fPid = p->PdgCode();
+ newp->fLab = nents;
+ Int_t moi = p->GetMother();
+ if (moi>=0) {
+ const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
+ moi = mmc->GetUniqueID();
+ }
+ newp->fMo = moi;
+ p->SetUniqueID(nents);
+
+ Int_t nref = p->GetNumberOfTrackReferences();
+ if (nref>0) {
+ AliTrackReference *ref = p->GetTrackReference(nref-1);
+ if (ref) {
+ newp->fDet = ref->DetectorId();
+ }
+ }
+
+ if (d1<0)
+ return;
+ if (d2<0)
+ d2=d1;
+ for (Int_t i=d1;i<=d2;++i) {
+ AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
+ if (dmc->P()<0.01)
+ continue;
+ ProcessDaughters(dmc,i,arr);
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff