#include "AliMCParticle.h"
#include "AliMCEvent.h"
#include "AliAnalysisHelperJetTasks.h"
-#include "AliAODMCParticle.h"
-
-using std::cout;
-using std::endl;
+#include "AliAODTrack.h"
ClassImp(AliMuonEffMC)
//________________________________________________________________________
AliMuonEffMC::AliMuonEffMC() :
- AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
- fHEventStat(0), fHXsec(0), fHTrials(0), fHEvt(0x0), fIsMc(kTRUE), fIsPythia(kFALSE), fMDProcess(kFALSE), fFeynmanX(kFALSE), fScatFX(kFALSE), fZvProcess(kTRUE),
- fIsCutStudy(kFALSE), fIsFPM(kTRUE), fCentralityEstimator("V0M"), fNEtaBins(15), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12), fNPBins(150), fChiSquareNormCut(5.0),
- fHMuonParGen(0x0), fHMuonParGenSec(0x0), fHMuonParGenFPM(0x0), fHMuonParGenSecFPM(0x0), fHMuonParGenP(0x0), fHMuonDetGenP(0x0), fHMuonDetRecP(0x0),
- fHFXu(0), fHFXantiu(0), fHFXd(0), fHFXantid(0), fHFXg(0), fHFXetc(0),
- fHaFx(0), fHbFx(0), fHabFxRatio(0), fHabDeltaFx(0), fHabRelDeltaFx(0),
- fHZvRvPrim(0x0), fHZvRvSec(0x0), fHXvYvPrim(0x0), fHXvYvSec(0x0)
+ AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99),
+ fOutputList(0x0), fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fPlotMode(0), fMuonCutMask(0), fMuonTrackCuts(0x0),
+ fCentralityEstimator("V0M"), fNEtaBins(100), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12),
+ fHFM(0x0), fHPP(0x0), fHMuFM(0x0), fHMuFMrec(0x0), fHMuPP(0x0), fHMuPPrec(0x0), fHMuRec(0x0)
{
// Constructor
//DefineInput(0, TChain::Class());
//DefineOutput(1, TList::Class());
- for(Int_t i=0; i<6; i++)
- {
- fHMuonDetGen[i] = NULL;
- fHMuonDetRec[i] = NULL;
- fHHadDetRec[i] = NULL;
- fHSecDetRec[i] = NULL;
- }
- for(Int_t i=0; i<4; i++)
- {
- fHMuonParGenV[i] = NULL;
- fHMuonParGenSecV[i] = NULL;
- fHMuonDetRecV[i] = NULL;
- fHMuMotherGenPt[i] = NULL;
- fHMuMotherRecPt[i] = NULL;
- fHMuMotherGenPhi[i] = NULL;
- fHMuMotherRecPhi[i] = NULL;
- fHMuMotherGenEta[i] = NULL;
- fHMuMotherRecEta[i] = NULL;
- fHMuDCA[i] = NULL;
- fHMuMotherGenPtSec[i] = NULL;
- fHMuMotherRecPtSec[i] = NULL;
- fHMuMotherGenPhiSec[i] = NULL;
- fHMuMotherRecPhiSec[i] = NULL;
- fHMuMotherGenEtaSec[i] = NULL;
- fHMuMotherRecEtaSec[i] = NULL;
- fHMuDCASec[i] = NULL;
- fHMuZv[i] = NULL;
- fHMuRelZv[i] = NULL;
- for(Int_t j=0; j<3; j++)
- {
- fHMuMohterPhiDifGen[i][j] = NULL;
- fHMuMohterPhiDifRec[i][j] = NULL;
- fHMuMohterEtaDifGen[i][j] = NULL;
- fHMuMohterEtaDifRec[i][j] = NULL;
- fHMuMohterPhiDifGenSec[i][j] = NULL;
- fHMuMohterPhiDifRecSec[i][j] = NULL;
- fHMuMohterEtaDifGenSec[i][j] = NULL;
- fHMuMohterEtaDifRecSec[i][j] = NULL;
- }
- }
- for(Int_t i=0; i<2; i++)
- {
- for(Int_t j=0; j<3; j++)
- {
- fHFXmuonP[i][j] = NULL;
- fHFXmuonM[i][j] = NULL;
- }
- }
- for(Int_t i=0; i<3; i++)
- {
- fHabFxMu[i] = NULL;
- fHabFxRatioMu[i] = NULL;
- fHabDeltaFxMu[i] = NULL;
- fHabRelDeltaFxMu[i] = NULL;
- }
}
//________________________________________________________________________
AliMuonEffMC::AliMuonEffMC(const char *name) :
- AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
- fHEventStat(0), fHXsec(0), fHTrials(0), fHEvt(0x0), fIsMc(kTRUE), fIsPythia(kFALSE), fMDProcess(kFALSE), fFeynmanX(kFALSE), fScatFX(kFALSE), fZvProcess(kTRUE),
- fIsCutStudy(kFALSE), fIsFPM(kTRUE), fCentralityEstimator("V0M"), fNEtaBins(15), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12), fNPBins(150), fChiSquareNormCut(5.0),
- fHMuonParGen(0x0), fHMuonParGenSec(0x0), fHMuonParGenFPM(0x0), fHMuonParGenSecFPM(0x0), fHMuonParGenP(0x0), fHMuonDetGenP(0x0), fHMuonDetRecP(0x0),
- fHFXu(0), fHFXantiu(0), fHFXd(0), fHFXantid(0), fHFXg(0), fHFXetc(0),
- fHaFx(0), fHbFx(0), fHabFxRatio(0), fHabDeltaFx(0), fHabRelDeltaFx(0),
- fHZvRvPrim(0x0), fHZvRvSec(0x0), fHXvYvPrim(0x0), fHXvYvSec(0x0)
+ AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99),
+ fOutputList(0x0), fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fPlotMode(0), fMuonCutMask(0), fMuonTrackCuts(0x0),
+ fCentralityEstimator("V0M"), fNEtaBins(100), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12),
+ fHFM(0x0), fHPP(0x0), fHMuFM(0x0), fHMuFMrec(0x0), fHMuPP(0x0), fHMuPPrec(0x0), fHMuRec(0x0)
{
// Constructor
- for(Int_t i=0; i<6; i++)
- {
- fHMuonDetGen[i] = NULL;
- fHMuonDetRec[i] = NULL;
- fHHadDetRec[i] = NULL;
- fHSecDetRec[i] = NULL;
- }
- for(Int_t i=0; i<4; i++)
- {
- fHMuonParGenV[i] = NULL;
- fHMuonParGenSecV[i] = NULL;
- fHMuonDetRecV[i] = NULL;
- fHMuMotherGenPt[i] = NULL;
- fHMuMotherRecPt[i] = NULL;
- fHMuMotherGenPhi[i] = NULL;
- fHMuMotherRecPhi[i] = NULL;
- fHMuMotherGenEta[i] = NULL;
- fHMuMotherRecEta[i] = NULL;
- fHMuDCA[i] = NULL;
- fHMuMotherGenPtSec[i] = NULL;
- fHMuMotherRecPtSec[i] = NULL;
- fHMuMotherGenPhiSec[i] = NULL;
- fHMuMotherRecPhiSec[i] = NULL;
- fHMuMotherGenEtaSec[i] = NULL;
- fHMuMotherRecEtaSec[i] = NULL;
- fHMuDCASec[i] = NULL;
- fHMuZv[i] = NULL;
- fHMuRelZv[i] = NULL;
- for(Int_t j=0; j<3; j++)
- {
- fHMuMohterPhiDifGen[i][j] = NULL;
- fHMuMohterPhiDifRec[i][j] = NULL;
- fHMuMohterEtaDifGen[i][j] = NULL;
- fHMuMohterEtaDifRec[i][j] = NULL;
- fHMuMohterPhiDifGenSec[i][j] = NULL;
- fHMuMohterPhiDifRecSec[i][j] = NULL;
- fHMuMohterEtaDifGenSec[i][j] = NULL;
- fHMuMohterEtaDifRecSec[i][j] = NULL;
- }
- }
- for(Int_t i=0; i<2; i++)
- {
- for(Int_t j=0; j<3; j++)
- {
- fHFXmuonP[i][j] = NULL;
- fHFXmuonM[i][j] = NULL;
- }
- }
- for(Int_t i=0; i<3; i++)
- {
- fHabFxMu[i] = NULL;
- fHabFxRatioMu[i] = NULL;
- fHabDeltaFxMu[i] = NULL;
- fHabRelDeltaFxMu[i] = NULL;
- }
DefineInput(0, TChain::Class());
DefineOutput(1, TList::Class());
}
fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt"); //!Global Trigger Single plus Low p_T
fOutputList->Add(fHEventStat);
- if(fIsPythia)
- {
- fHXsec = new TH1F("fHXsec", "Cross section from pyxsec.root", 1, 0, 1);
- fHXsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
- fOutputList->Add(fHXsec);
-
- fHTrials = new TH1F("fHTrials", "Number of Trials", 1, 0, 1);
- fHTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
- fOutputList->Add(fHTrials);
- }
-
fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
fOutputList->Add(fHEvt);
+ fMuonTrackCuts = new AliMuonTrackCuts("StdMuonCuts","StdMuonCuts");
+ fMuonTrackCuts->SetCustomParamFromRun(197388,"muon_pass2"); // for LHC13b,c,d,e,f
+ fMuonTrackCuts->SetFilterMask(fMuonCutMask);
+ AliInfo(Form(" using muon track cuts with bit %u\n", fMuonCutMask));
+
+ // Define THn's
Int_t iTrackBin[6];
Double_t* trackBins[6];
const char* trackAxisTitle[6];
- Int_t iTrackBinP[6];
- Double_t* trackBinsP[6];
- const char* trackAxisTitleP[6];
-
// eta
Double_t etaBins[fNEtaBins+1];
- for(Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
+ for(Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-8.0 + 16.0/fNEtaBins*i); }
iTrackBin[0] = fNEtaBins;
- iTrackBinP[0] = fNEtaBins;
trackBins[0] = etaBins;
- trackBinsP[0] = etaBins;
trackAxisTitle[0] = "#eta";
- trackAxisTitleP[0] = "#eta";
- // p_T
+ // p_T
Double_t pTBins[fNpTBins+1];
for(Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(5.0/fNpTBins * i); }
iTrackBin[1] = fNpTBins;
trackBins[1] = pTBins;
trackAxisTitle[1] = "p_{T} (GeV/c)";
- // P
- Double_t PBins[fNPBins+1];
- for(Int_t i=0; i<=fNPBins; i++) { PBins[i] = (Double_t)(150.0/fNPBins * i); }
- iTrackBinP[1] = fNPBins;
- trackBinsP[1] = PBins;
- trackAxisTitleP[1] = "P (GeV/c)";
-
- // centrality
+ // centrality/multiplicity
Double_t CentBins[fNCentBins+1];
for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins * i); }
iTrackBin[2] = fNCentBins;
- iTrackBinP[2] = fNCentBins;
trackBins[2] = CentBins;
- trackBinsP[2] = CentBins;
- trackAxisTitle[2] = "Cent";
- trackAxisTitleP[2] = "Cent";
-
- // Z-vertex
- Double_t ZvtxBins[fNZvtxBins+1];
- for(Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins * i); }
- iTrackBin[3] = fNZvtxBins;
- iTrackBinP[3] = fNZvtxBins;
- trackBins[3] = ZvtxBins;
- trackBinsP[3] = ZvtxBins;
- trackAxisTitle[3] = "Zvtx";
- trackAxisTitleP[3] = "Zvtx";
+ trackAxisTitle[2] = "Mult";
// phi
Double_t phiBins[fNPhiBins+1];
for(Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
- iTrackBin[4] = fNPhiBins;
- iTrackBinP[4] = fNPhiBins;
- trackBins[4] = phiBins;
- trackBinsP[4] = phiBins;
- trackAxisTitle[4] = "#phi";
- trackAxisTitleP[4] = "#phi";
-
- // charge
- Double_t chargeBins[3] = {-10.0, 0, 10.0};
- iTrackBin[5] = 2;
- iTrackBinP[5] = 2;
- trackBins[5] = chargeBins;
- trackBinsP[5] = chargeBins;
- trackAxisTitle[5] = "charge";
- trackAxisTitleP[5] = "charge";
+ iTrackBin[3] = fNPhiBins;
+ trackBins[3] = phiBins;
+ trackAxisTitle[3] = "#phi";
+
+ // charge
+ Double_t chargeBins[4] = {-10.0, -0.5, 0.5, 10.0};
+ iTrackBin[4] = 3;
+ trackBins[4] = chargeBins;
+ trackAxisTitle[4] = "charge";
+
+ // species
+ Double_t Species[21];
+ for(Int_t iSpe=0; iSpe<=20; iSpe++) { Species[iSpe] = (Double_t)iSpe; }
+ iTrackBin[5] = 20;
+ trackBins[5] = Species;
+ trackAxisTitle[5] = "Species";
- const char *cutlabel[6] = {"NoCut", "Muon", "Trigger", "eta", "ThetaAbs", "ChiS"};
if(fIsMc)
{
// THn for tracking efficiency
- fHMuonParGen = new THnF("fHMuonParGen", "", 6, iTrackBin, 0, 0);
- for (Int_t i=0; i<6; i++)
+ if(fPlotMode==0)
{
- fHMuonParGen->SetBinEdges(i, trackBins[i]);
- fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
- }
- fHMuonParGen->Sumw2();
- fOutputList->Add(fHMuonParGen);
-
- fHMuonParGenSec = (THnF*) fHMuonParGen->Clone("fHMuonParGenSec");
- fHMuonParGenSec->Sumw2();
- fOutputList->Add(fHMuonParGenSec);
- if(fIsFPM)
- {
- fHMuonParGenFPM = (THnF*) fHMuonParGen->Clone("fHMuonParGenFPM");
- fHMuonParGenFPM->Sumw2();
- fOutputList->Add(fHMuonParGenFPM);
-
- fHMuonParGenSecFPM = (THnF*) fHMuonParGen->Clone("fHMuonParGenSecFPM");
- fHMuonParGenSecFPM->Sumw2();
- fOutputList->Add(fHMuonParGenSecFPM);
- }
- if(fIsCutStudy)
- {
- for(Int_t i=0; i<6; i++)
+ fHFM = new THnF("fHFM", "", 6, iTrackBin, 0, 0);
+ for (Int_t i=0; i<6; i++)
{
- fHMuonDetGen[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonDetGen_%s",cutlabel[i]));
- fHMuonDetGen[i]->Sumw2();
- fOutputList->Add(fHMuonDetGen[i]);
-
- fHMuonDetRec[i] = (THnF*) fHMuonParGen->Clone(Form("fHMuonDetRec_%s",cutlabel[i]));
- fHMuonDetRec[i]->Sumw2();
- fOutputList->Add(fHMuonDetRec[i]);
-
- fHHadDetRec[i] = (THnF*) fHMuonParGen->Clone(Form("fHHadDetRec_%s",cutlabel[i]));
- fHHadDetRec[i]->Sumw2();
- fOutputList->Add(fHHadDetRec[i]);
-
- fHSecDetRec[i] = (THnF*) fHMuonParGen->Clone(Form("fHSecDetRec_%s",cutlabel[i]));
- fHSecDetRec[i]->Sumw2();
- fOutputList->Add(fHSecDetRec[i]);
+ fHFM->SetBinEdges(i, trackBins[i]);
+ fHFM->GetAxis(i)->SetTitle(trackAxisTitle[i]);
}
+ fOutputList->Add(fHFM);
}
- else
- {
- fHMuonDetGen[0] = (THnF*)fHMuonParGen->Clone(Form("fHMuonDetGen"));
- fHMuonDetGen[0]->Sumw2();
- fOutputList->Add(fHMuonDetGen[0]);
-
- fHMuonDetRec[0] = (THnF*) fHMuonParGen->Clone(Form("fHMuonDetRec"));
- fHMuonDetRec[0]->Sumw2();
- fOutputList->Add(fHMuonDetRec[0]);
-
- fHHadDetRec[0] = (THnF*) fHMuonParGen->Clone(Form("fHHadDetRec"));
- fHHadDetRec[0]->Sumw2();
- fOutputList->Add(fHHadDetRec[0]);
-
- fHSecDetRec[0] = (THnF*) fHMuonParGen->Clone(Form("fHSecDetRec"));
- fHSecDetRec[0]->Sumw2();
- fOutputList->Add(fHSecDetRec[0]);
- }
-
- const char *vertexlabel[4] = {"0_10", "10_90", "90_503", "503"};
- for(Int_t i=0; i<4; i++)
- {
- fHMuonParGenV[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonParGenV_%s",vertexlabel[i]));
- fHMuonParGenV[i]->Sumw2();
- fOutputList->Add(fHMuonParGenV[i]);
-
- fHMuonParGenSecV[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonParGenSecV_%s",vertexlabel[i]));
- fHMuonParGenSecV[i]->Sumw2();
- fOutputList->Add(fHMuonParGenSecV[i]);
-
- fHMuonDetRecV[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonDetRecV_%s",vertexlabel[i]));
- fHMuonDetRecV[i]->Sumw2();
- fOutputList->Add(fHMuonDetRecV[i]);
- }
-
- fHMuonParGenP = new THnF("fHMuonParGenP", "", 6, iTrackBinP, 0, 0);
- for (Int_t i=0; i<6; i++)
+ else if(fPlotMode==1)
{
- fHMuonParGenP->SetBinEdges(i, trackBinsP[i]);
- fHMuonParGenP->GetAxis(i)->SetTitle(trackAxisTitleP[i]);
- }
- fHMuonParGenP->Sumw2();
- fOutputList->Add(fHMuonParGenP);
-
- fHMuonDetGenP = (THnF*) fHMuonParGenP->Clone("fHMuonDetGenP");
- fHMuonDetGenP->Sumw2();
- fOutputList->Add(fHMuonDetGenP);
-
- fHMuonDetRecP = (THnF*) fHMuonParGenP->Clone("fHMuonDetRecP");
- fHMuonDetRecP->Sumw2();
- fOutputList->Add(fHMuonDetRecP);
-
- const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
- const char *MuPt[3] = {"0005","0520","2040"};
- if(fZvProcess)
- {
- for(Int_t i=0; i<4; i++)
+ fHPP = new THnF("fHPP", "", 6, iTrackBin, 0, 0);
+ for (Int_t i=0; i<6; i++)
{
- fHMuZv[i] = new TH1F(Form("fHMuZv_%s",MotherSpecies[i]), ";Z_{V}(cm)", 600, -500, 100);
- fOutputList->Add(fHMuZv[i]);
- fHMuRelZv[i] = new TH1F(Form("fHMuRelZv_%s",MotherSpecies[i]), ";|Z_{V,muon}-Zvtx|(cm)", 500, 0, 500);
- fOutputList->Add(fHMuRelZv[i]);
+ fHPP->SetBinEdges(i, trackBins[i]);
+ fHPP->GetAxis(i)->SetTitle(trackAxisTitle[i]);
}
+ fOutputList->Add(fHPP);
}
- if(fMDProcess)
+
+ else if(fPlotMode==2)
{
- for(Int_t i=0; i<4; i++)
+ fHMuFM = new THnF("fHMuFM", "", 6, iTrackBin, 0, 0);
+ for (Int_t i=0; i<6; i++)
{
- fHMuMotherGenPt[i] = new TH2F(Form("fHMuMotherGenPt_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
- fOutputList->Add(fHMuMotherGenPt[i]);
- fHMuMotherRecPt[i] = new TH2F(Form("fHMuMotherRecPt_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
- fOutputList->Add(fHMuMotherRecPt[i]);
- fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
- fOutputList->Add(fHMuMotherGenPhi[i]);
- fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
- fOutputList->Add(fHMuMotherRecPhi[i]);
- fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
- fOutputList->Add(fHMuMotherGenEta[i]);
- fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
- fOutputList->Add(fHMuMotherRecEta[i]);
- fHMuDCA[i] = new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
- fOutputList->Add(fHMuDCA[i]);
-
- fHMuMotherGenPtSec[i] = new TH2F(Form("fHMuMotherGenPtSec_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
- fOutputList->Add(fHMuMotherGenPtSec[i]);
- fHMuMotherRecPtSec[i] = new TH2F(Form("fHMuMotherRecPtSec_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
- fOutputList->Add(fHMuMotherRecPtSec[i]);
- fHMuMotherGenPhiSec[i] = new TH2F(Form("fHMuMotherGenPhiSec_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
- fOutputList->Add(fHMuMotherGenPhiSec[i]);
- fHMuMotherRecPhiSec[i] = new TH2F(Form("fHMuMotherRecPhiSec_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
- fOutputList->Add(fHMuMotherRecPhiSec[i]);
- fHMuMotherGenEtaSec[i] = new TH2F(Form("fHMuMotherGenEtaSec_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
- fOutputList->Add(fHMuMotherGenEtaSec[i]);
- fHMuMotherRecEtaSec[i] = new TH2F(Form("fHMuMotherRecEtaSec_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
- fOutputList->Add(fHMuMotherRecEtaSec[i]);
- fHMuDCASec[i] = new TH1F(Form("fHMuDCASec_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
- fOutputList->Add(fHMuDCASec[i]);
-
- for(Int_t j=0; j<3; j++)
- {
- fHMuMohterPhiDifGen[i][j] = new TH1F(Form("fHMuMohterPhiDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
- fOutputList->Add(fHMuMohterPhiDifGen[i][j]);
- fHMuMohterPhiDifRec[i][j] = new TH1F(Form("fHMuMohterPhiDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
- fOutputList->Add(fHMuMohterPhiDifRec[i][j]);
- fHMuMohterEtaDifGen[i][j] = new TH1F(Form("fHMuMohterEtaDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
- fOutputList->Add(fHMuMohterEtaDifGen[i][j]);
- fHMuMohterEtaDifRec[i][j] = new TH1F(Form("fHMuMohterEtaDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
- fOutputList->Add(fHMuMohterEtaDifRec[i][j]);
-
- fHMuMohterPhiDifGenSec[i][j] = new TH1F(Form("fHMuMohterPhiDifGenSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
- fOutputList->Add(fHMuMohterPhiDifGenSec[i][j]);
- fHMuMohterPhiDifRecSec[i][j] = new TH1F(Form("fHMuMohterPhiDifRecSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
- fOutputList->Add(fHMuMohterPhiDifRecSec[i][j]);
- fHMuMohterEtaDifGenSec[i][j] = new TH1F(Form("fHMuMohterEtaDifGenSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
- fOutputList->Add(fHMuMohterEtaDifGenSec[i][j]);
- fHMuMohterEtaDifRecSec[i][j] = new TH1F(Form("fHMuMohterEtaDifRecSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
- fOutputList->Add(fHMuMohterEtaDifRecSec[i][j]);
- }
+ fHMuFM->SetBinEdges(i, trackBins[i]);
+ fHMuFM->GetAxis(i)->SetTitle(trackAxisTitle[i]);
}
- fHZvRvPrim = new TH2F("fHZvRvPrim", "", 300, -500, 100, 200, 0, 800);
- fOutputList->Add(fHZvRvPrim);
-
- fHZvRvSec = new TH2F("fHZvRvSec", "", 300, -500, 100, 200, 0, 800);
- fOutputList->Add(fHZvRvSec);
+ fOutputList->Add(fHMuFM);
- fHXvYvPrim = new TH2F("fHXvYvPrim", "", 200, -500, 500, 200, -500, 500);
- fOutputList->Add(fHXvYvPrim);
-
- fHXvYvSec = new TH2F("fHXvYvSec", "", 200, -500, 500, 200, -500, 500);
- fOutputList->Add(fHXvYvSec);
- }
- if(fFeynmanX)
- {
- fHFXu = new TH1F("fHFXu",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXu);
-
- fHFXantiu = new TH1F("fHFXantiu",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXantiu);
-
- fHFXd = new TH1F("fHFXd",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXd);
-
- fHFXantid = new TH1F("fHFXantid",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXantid);
-
- fHFXg = new TH1F("fHFXg",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXg);
-
- fHFXetc = new TH1F("fHFXetc",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXetc);
-
- const char* muonptbin[2] = {"02","2"};
- const char *muonetabin[3] = {"4035","3530","3025"};
- for(Int_t i=0; i<2; i++)
- {
- for(Int_t j=0; j<3; j++)
- {
- fHFXmuonP[i][j] = new TH1F(Form("fHFXmuonP_%s_%s",muonptbin[i], muonetabin[j]),";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXmuonP[i][j]);
- fHFXmuonM[i][j] = new TH1F(Form("fHFXmuonM_%s_%s",muonptbin[i], muonetabin[j]),";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHFXmuonM[i][j]);
- }
- }
+ fHMuFMrec = (THnF*)fHMuFM->Clone("fHMuFMrec");
+ fOutputList->Add(fHMuFMrec);
}
- if(fScatFX)
+ else if(fPlotMode==3)
{
- fHaFx = new TH1F("fHaFx",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHaFx);
- fHbFx = new TH1F("fHbFx",";x_{F}",1000, 0.0, 1.0);
- fOutputList->Add(fHbFx);
- fHabFxRatio = new TH1F("fHabFxRatio",";x_{F,a}/x_{F,b}",1000, 0.0, 5.0);
- fOutputList->Add(fHabFxRatio);
- fHabDeltaFx = new TH1F("fHabDeltaFx",";#Deltax_{F}",1000, -2.0, 2.0);
- fOutputList->Add(fHabDeltaFx);
- fHabRelDeltaFx = new TH1F("fHabRelDeltaFx",";(x_{a}-x_{b})/(x_{a}+x{b})",1000, -1.0, 1.0);
- fOutputList->Add(fHabRelDeltaFx);
-
- const char* MuonPt[3] = {"Mu01","Mu12", "Mu2"};
- for(Int_t i=0; i<3; i++)
+ fHMuPP = new THnF("fHMuPP", "", 6, iTrackBin, 0, 0);
+ for (Int_t i=0; i<6; i++)
{
- Int_t NabBin = 22;
- Double_t abBins[22+1] = { 0.0, 0.000001, 0.000002, 0.000004, 0.00001, 0.00002, 0.00004, 0.0001, 0.0002, 0.0004, 0.001, 0.002, 0.004, 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 1};
-
- fHabFxMu[i] = new TH2F(Form("fHabFxMu_%s",MuonPt[i]),";x_{F,a};x_{F,b}",NabBin, abBins, NabBin, abBins);
- fOutputList->Add(fHabFxMu[i]);
-
- fHabFxRatioMu[i] = new TH1F(Form("fHabFxRatioMu_%s",MuonPt[i]),";x_{F,a}/x_{F,b}",1000, 0.0, 5.0);
- fOutputList->Add(fHabFxRatioMu[i]);
-
- fHabDeltaFxMu[i] = new TH1F(Form("fHabDeltaFxMu_%s",MuonPt[i]),";#Deltax_{F}",1000, -2.0, 2.0);
- fOutputList->Add(fHabDeltaFxMu[i]);
-
- fHabRelDeltaFxMu[i] = new TH1F(Form("fHabRelDeltaFxMu_%s",MuonPt[i]),";(x_{a}-x_{b})/(x_{a}+x{b})",1000, -1.0, 1.0);
- fOutputList->Add(fHabRelDeltaFxMu[i]);
+ fHMuPP->SetBinEdges(i, trackBins[i]);
+ fHMuPP->GetAxis(i)->SetTitle(trackAxisTitle[i]);
}
+ fOutputList->Add(fHMuPP);
+
+ fHMuPPrec = (THnF*)fHMuPP->Clone("fHMuPPrec");
+ fOutputList->Add(fHMuPPrec);
}
}
else
{
- fHMuonDetRec[0] = new THnF(Form("fHMuonDetRec_%s", cutlabel[0]), "", 6, iTrackBin, 0, 0);
- for (Int_t i=0; i<6; i++)
- {
- fHMuonDetRec[0]->SetBinEdges(i, trackBins[i]);
- fHMuonDetRec[0]->GetAxis(i)->SetTitle(trackAxisTitle[i]);
- }
- fHMuonDetRec[0]->Sumw2();
- fOutputList->Add(fHMuonDetRec[0]);
-
- for(Int_t i=1; i<6; i++)
+ if(fPlotMode==4)
{
-
- fHMuonDetRec[i] = (THnF*) fHMuonDetRec[0]->Clone(Form("fHMuonDetRec_%s",cutlabel[i]));
- fHMuonDetRec[i]->Sumw2();
- fOutputList->Add(fHMuonDetRec[i]);
- }
-
- fHMuonDetRecP = new THnF("fHMuonDetRecP", "", 6, iTrackBinP, 0, 0);
- for (Int_t i=0; i<6; i++)
- {
- fHMuonDetRecP->SetBinEdges(i, trackBinsP[i]);
- fHMuonDetRecP->GetAxis(i)->SetTitle(trackAxisTitleP[i]);
+ fHMuRec = new THnF("fHMuRec", "", 6, iTrackBin, 0, 0);
+ for (Int_t i=0; i<6; i++)
+ {
+ fHMuRec->SetBinEdges(i, trackBins[i]);
+ fHMuRec->GetAxis(i)->SetTitle(trackAxisTitle[i]);
+ }
+ fOutputList->Add(fHMuRec);
}
- fHMuonDetRecP->Sumw2();
- fOutputList->Add(fHMuonDetRecP);
}
PostData(1, fOutputList);
}
-//________________________________________________________________________
-Bool_t AliMuonEffMC::Notify()
-{
- // Implemented Notify() to read the cross sections
- // and number of trials from pyxsec.root
- if(fIsPythia)
- {
- fHEventStat->Fill(2.5);
-
- TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
- Float_t xsection = 0;
- Float_t ftrials = 1;
-
- if(tree){
- TFile *curfile = tree->GetCurrentFile();
- if (!curfile) {
- Error("Notify","No current file");
- return kFALSE;
- }
-
- AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
- fHXsec->Fill("<#sigma>",xsection);
- fHTrials->Fill("#sum{ntrials}",ftrials);
- }
- }
- return kTRUE;
-}
//________________________________________________________________________
void AliMuonEffMC::UserExec(Option_t *)
{
{
fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
- ntrks = fAOD->GetNTracks();
+ ntrks = fAOD->GetNumberOfTracks();
}
else
{
if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
}
- if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) { //AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
- return; }
- if (fCentrality > 100. || fCentrality < -1.5) { //AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
- return; }
+ if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) {
+ //AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
+ return;
+ }
+ if (fCentrality > 100. || fCentrality < -1.5) {
+ //AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
+ return;
+ }
if(fCentrality < 0) fCentrality = 0.5; //ad hoc centrality for pp
// Fill Event histogram
if (trigword & 0x1000) fHEventStat->Fill(14.5);
if (trigword & 0x2000) fHEventStat->Fill(15.5);
if (trigword & 0x4000) fHEventStat->Fill(16.5);
-
- if(fIsMc)
+
+ if(fIsMc && (fPlotMode==0 || fPlotMode==1) && fESD)
{
- // generated level loop
- for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
+ Double_t MEta = 0.0;
+ Double_t MPt = 0.0;
+ Double_t MPhi = 0.0;
+ Double_t MCharge = 0.0;
+ Double_t MSpecies = 0.0;
+
+ for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++) // generated level loop
{
- Int_t FirstPrimMother = 0;
- Double_t FirstPrimZvtx = 0.0;
- Int_t Zvtxbin = 0;
- if(fAOD)
+ AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(ipart);
+ if(fPlotMode==0)
{
- AliAODMCParticle *AodMcParticle = (AliAODMCParticle*)fMC->GetTrack(ipart);
- if(AodMcParticle->Eta() < -4.0 || AodMcParticle->Eta() > -2.5) continue;
-
- if(TMath::Abs(AodMcParticle->PdgCode())==13)
- {
- FirstPrimMother = GetFirstPrimaryMother(ipart);
- FirstPrimZvtx = GetFirstPrimaryVertex(ipart);
- Zvtxbin = GetZVertexBin(FirstPrimZvtx);
-
- Double_t fillArrayParGen[6] = { AodMcParticle->Eta(), AodMcParticle->Pt(), fCentrality, fZVertex, AodMcParticle->Phi(), AodMcParticle->Charge() };
- if(AodMcParticle->IsPrimary())
- {
- fHMuonParGen->Fill(fillArrayParGen);
- fHMuonParGenV[Zvtxbin]->Fill(fillArrayParGen);
- }
- else
- {
- fHMuonParGenSec->Fill(fillArrayParGen);
- fHMuonParGenSecV[Zvtxbin]->Fill(fillArrayParGen);
- }
- Double_t fillArrayParGenP[6] = { AodMcParticle->Eta(), AodMcParticle->P(), fCentrality, fZVertex, AodMcParticle->Phi(), AodMcParticle->Charge() };
- if(AodMcParticle->IsPrimary()) fHMuonParGenP->Fill(fillArrayParGenP);
- }
+ if(ipart < 211 || McParticle->GetMother()!=-1) continue;
+ MEta = McParticle->Eta();
+ MPt = McParticle->Pt();
+ MPhi = McParticle->Phi();
+ MCharge = McParticle->Charge();
+ MSpecies = GetSpecies(McParticle->PdgCode());
}
- else if(fESD)
- {
- AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(ipart);
- if(McParticle->Eta() < -4.0 || McParticle->Eta() > -2.5) continue;
-
- if(TMath::Abs(McParticle->PdgCode())==13)
- {
- FirstPrimMother = GetFirstPrimaryMother(ipart);
- FirstPrimZvtx = GetFirstPrimaryVertex(ipart);
- Zvtxbin = GetZVertexBin(FirstPrimZvtx);
-
- Double_t fillArrayParGen[6] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi(), McParticle->Charge() };
- AliMCParticle *FPrimaryMother = (AliMCParticle*)fMC->GetTrack(FirstPrimMother);
- Double_t fillArrayParGenFPM[6] = { FPrimaryMother->Eta(), FPrimaryMother->Pt(), fCentrality, fZVertex, FPrimaryMother->Phi(), FPrimaryMother->Charge() };
-
- if(McParticle->GetMother()< fStack->GetNprimary())
- {
- fHMuonParGen->Fill(fillArrayParGen);
- fHMuonParGenV[Zvtxbin]->Fill(fillArrayParGen);
- if(fIsFPM) fHMuonParGenFPM->Fill(fillArrayParGenFPM);
- }
- else
- {
- fHMuonParGenSec->Fill(fillArrayParGen);
- fHMuonParGenSecV[Zvtxbin]->Fill(fillArrayParGen);
- if(fIsFPM) fHMuonParGenSecFPM->Fill(fillArrayParGenFPM);
- }
- Double_t fillArrayParGenP[6] = { McParticle->Eta(), McParticle->P(), fCentrality, fZVertex, McParticle->Phi(), McParticle->Charge() };
- if(McParticle->GetMother()< fStack->GetNprimary()) fHMuonParGenP->Fill(fillArrayParGenP);
- }
+ else if(fPlotMode==1)
+ {
+ if(!fMC->IsPhysicalPrimary(ipart)) continue;
+ MEta = McParticle->Eta();
+ MPt = McParticle->Pt();
+ MPhi = McParticle->Phi();
+ MCharge = McParticle->Charge();
+ MSpecies = GetSpecies(McParticle->PdgCode());
}
- }
- }
-
- // reconstructed level loop
- for (Int_t iTrack = 0; iTrack<ntrks; iTrack++)
+ Double_t fillArrayM[6] = { MEta, MPt, fCentrality, MPhi, MCharge, MSpecies };
+ if(fPlotMode==0) fHFM->Fill(fillArrayM);
+ else if(fPlotMode==1) fHPP->Fill(fillArrayM);
+ }// end of generated level loop
+ }
+
+ if((fPlotMode==2 || fPlotMode==3 || fPlotMode==4))
{
- Int_t label = 0;
- // reconstructed track variables
- Double_t trackpt = 0;
- Double_t tracketa = 0;
- Double_t trackphi = 0;
- Double_t trackcharge = 0;
- Double_t trackp = 0;
- Double_t dcavalue = 0;
- Int_t cutNum = 0;
- Bool_t isprimary = kFALSE;
- // reconstructed track's matched truth particle variables
- Double_t mcpt = 0;
- Double_t mceta = 0;
- Double_t mcphi = 0;
- Double_t mcp = 0;
- Double_t mccharge = 0;
- Double_t mcZv = 0;
- // first primary mother variables
- Int_t motherlabel =0;
- Int_t motherpdg = 0;
- Double_t motherpt = 0;
- Double_t mothereta = 0;
- Double_t motherphi = 0;
- Double_t motherXv = 0;
- Double_t motherYv = 0;
- Double_t motherZv = 0;
-
- if(fAOD)
+ for (Int_t iTrack = 0; iTrack<ntrks; iTrack++) // reconstructed level loop
{
- AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
- if(muonTrack)
- {
- if(!(muonTrack->IsMuonTrack())) continue;
- cutNum = IsGoodMUONtrack(*muonTrack);
- trackpt = muonTrack->Pt();
- tracketa = muonTrack->Eta();
- trackphi = muonTrack->Phi();
- trackcharge = muonTrack->Charge();
- trackp = muonTrack->P();
- label = TMath::Abs(muonTrack->GetLabel());
- if (label>=fMC->GetNumberOfTracks()) {
- AliError(Form("Label %d larger than number of particles on stack %d\n",label,fMC->GetNumberOfTracks()));
- continue;
- }
- }
- }
- else if(fESD)
- {
- AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
- if(muonTrack)
- {
- cutNum = IsGoodMUONtrack(*muonTrack);
- trackpt = muonTrack->Pt();
- tracketa = muonTrack->Eta();
- trackphi = muonTrack->Phi();
- trackcharge = muonTrack->Charge();
- trackp = muonTrack->P();
- label = TMath::Abs(muonTrack->GetLabel());
- dcavalue = muonTrack->GetDCA();
- }
- }
- Double_t fillArrayDetRec[6] = { tracketa, trackpt, fCentrality, fZVertex, trackphi, trackcharge };
- Double_t fillArrayDetRecP[6] = { tracketa, trackp, fCentrality, fZVertex, trackphi, trackcharge };
- if(fIsCutStudy) { for(Int_t icut = 0; icut <= cutNum; icut++) { fHMuonDetRec[icut]->Fill(fillArrayDetRec); }}
- else { if(cutNum==5) fHMuonDetRec[0]->Fill(fillArrayDetRec); }
- if(cutNum==5) fHMuonDetRecP->Fill(fillArrayDetRecP);
-
- if(fIsMc)
- {
- if(fAOD)
- {
- AliAODMCParticle *aodMcParticle = (AliAODMCParticle*)fMC->GetTrack(label);
- isprimary = (aodMcParticle->IsPrimary()) ? kTRUE : kFALSE;
- if(TMath::Abs(aodMcParticle->PdgCode()) != 13)
- {
- if(fIsCutStudy){ for(Int_t icut=0; icut<=cutNum; icut++) fHHadDetRec[icut]->Fill(fillArrayDetRec); }
- else { if(cutNum==5) fHHadDetRec[0]->Fill(fillArrayDetRec); }
- }
- else if(!isprimary)
- {
- if(fIsCutStudy){ for(Int_t icut=0; icut<=cutNum; icut++) fHSecDetRec[icut]->Fill(fillArrayDetRec); }
- else { if(cutNum==5) fHSecDetRec[0]->Fill(fillArrayDetRec); }
- }
- mcpt = aodMcParticle->Pt();
- mceta = aodMcParticle->Eta();
- mcphi = aodMcParticle->Phi();
- mccharge = aodMcParticle->Charge();
- mcp = aodMcParticle->P();
- mcZv = aodMcParticle->Zv();
- motherlabel = aodMcParticle->GetMother();
- }
- else if(fESD)
+ Int_t label = 0;
+ Double_t MuEta = 0.0;
+ Double_t MuPt = 0.0;
+ Double_t MuPhi = 0.0;
+ Double_t MuCharge = 0.0;
+
+ Double_t MuMEta = 0.0;
+ Double_t MuMPt = 0.0;
+ Double_t MuMPhi = 0.0;
+ Double_t MuMCharge = 0.0;
+ Double_t MuMSpecies = 0.0;
+
+ if(fESD)
{
- AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(label);
- isprimary = (McParticle->GetMother()<fStack->GetNprimary()) ? kTRUE : kFALSE;
- if(TMath::Abs(McParticle->PdgCode())!=13)
- {
- if(fIsCutStudy) { for(Int_t icut=0; icut<=cutNum; icut++) fHHadDetRec[icut]->Fill(fillArrayDetRec); }
- else{ if(cutNum==5) fHHadDetRec[0]->Fill(fillArrayDetRec); }
- }
- if(!isprimary)
- {
- if(fIsCutStudy) { for(Int_t icut=0; icut<=cutNum; icut++) fHSecDetRec[icut]->Fill(fillArrayDetRec); }
- else{ if(cutNum==5) fHSecDetRec[0]->Fill(fillArrayDetRec); }
- }
- mcpt = McParticle->Pt();
- mceta = McParticle->Eta();
- mcphi = McParticle->Phi();
- mccharge = McParticle->Charge();
- mcp = McParticle->P();
- mcZv = McParticle->Zv();
-
- motherlabel = GetFirstPrimaryMother(label);
- Double_t primzvtx = GetFirstPrimaryVertex(label);
- Int_t reczvtxbin = GetZVertexBin(primzvtx);
- if(cutNum==5) fHMuonDetRecV[reczvtxbin]->Fill(fillArrayDetRec);
-
- if(motherlabel > -1)
+ AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
+ if(muonTrack)
{
- AliMCParticle *FPrimaryMother = (AliMCParticle*)fMC->GetTrack(motherlabel);
- motherpdg = TMath::Abs(FPrimaryMother->PdgCode());
- motherpt = FPrimaryMother->Pt();
- mothereta = FPrimaryMother->Eta();
- motherphi = FPrimaryMother->Phi();
- AliMCParticle *DaughtParticle = (AliMCParticle*)fMC->GetTrack(FPrimaryMother->GetFirstDaughter());
- motherXv = DaughtParticle->Xv();
- motherYv = DaughtParticle->Yv();
- motherZv = DaughtParticle->Zv();
+ if(!IsGoodMUONtrack(*muonTrack)) continue;
+ MuEta = muonTrack->Eta();
+ MuPt = muonTrack->Pt();
+ MuPhi = muonTrack->Phi();
+ MuCharge = muonTrack->Charge();
+ label = TMath::Abs(muonTrack->GetLabel());
+ if (label>=fMC->GetNumberOfTracks()) {
+ AliError(Form("Label %d larger than number of particles on stack %d\n",label,fMC->GetNumberOfTracks()));
+ continue;
+ }
}
}
- Double_t fillArrayDetGen[6] = { mceta, mcpt, fCentrality, fZVertex, mcphi, mccharge };
- if(fIsCutStudy){ for(Int_t icut=0; icut <= cutNum; icut++) fHMuonDetGen[icut]->Fill(fillArrayDetGen); }
- else {if(cutNum==5) fHMuonDetGen[0]->Fill(fillArrayDetGen); }
-
- Double_t fillArrayDetGenP[6] = { mceta, mcp, fCentrality, fZVertex, mcphi, mccharge };
- if(cutNum==5)fHMuonDetGenP->Fill(fillArrayDetGenP);
-
- Int_t motherbin = GetMotherBin(motherpdg);
- // muon Z-Vertex process
- if(fZvProcess && cutNum==5)
+ else if(fAOD)
{
- fHMuZv[motherbin]->Fill(mcZv);
- fHMuRelZv[motherbin]->Fill(TMath::Abs(mcZv-fZVertex));
+ AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
+ if(muonTrack)
+ {
+ if(!IsGoodMUONtrack(*muonTrack)) continue;
+ MuEta = muonTrack->Eta();
+ MuPt = muonTrack->Pt();
+ MuPhi = muonTrack->Phi();
+ MuCharge = muonTrack->Charge();
+ label = TMath::Abs(muonTrack->GetLabel());
+ if (label>=fMC->GetNumberOfTracks()) {
+ AliError(Form("Label %d larger than number of particles on stack %d\n",label,fMC->GetNumberOfTracks()));
+ continue;
+ }
+ }
}
- // mother-daughter kinematic relation
- if(fMDProcess && cutNum==5 && motherlabel>0)
+
+ if(fIsMc && fESD)
{
- if(isprimary)
+ if(fPlotMode==2)
{
- fHZvRvPrim->Fill(motherZv, TMath::Sqrt(motherXv*motherXv + motherYv*motherYv));
- fHXvYvPrim->Fill(motherXv, motherYv);
- }
- else
+ AliMCParticle* MuMparticle = (AliMCParticle*)fMC->GetTrack(GetFirstMother(label));
+ MuMEta = MuMparticle->Eta();
+ MuMPt = MuMparticle->Pt();
+ MuMPhi = MuMparticle->Phi();
+ MuMCharge = MuMparticle->Charge();
+ MuMSpecies = GetSpecies(MuMparticle->PdgCode());
+ }
+ else if(fPlotMode==3)
{
- fHZvRvSec->Fill(motherZv, TMath::Sqrt(motherXv*motherXv + motherYv*motherYv));
- fHXvYvSec->Fill(motherXv, motherYv);
+ AliMCParticle* MuMparticle = (AliMCParticle*)fMC->GetTrack(GetFirstPPMother(label));
+ MuMEta = MuMparticle->Eta();
+ MuMPt = MuMparticle->Pt();
+ MuMPhi = MuMparticle->Phi();
+ MuMCharge = MuMparticle->Charge();
+ MuMSpecies = GetSpecies(MuMparticle->PdgCode());
}
- MDProcess(isprimary, motherpdg, mcpt, mcphi, mceta, trackpt, trackphi, tracketa, motherpt, motherphi, mothereta, dcavalue);
- }
- }
- }
-
- if(fIsMc)
- {
- if(fFeynmanX) FeynmanX();
- if(fScatFX) ScatFX();
+ }// end of MC process
+ Double_t fillArrayMuRec[6] = { MuEta, MuPt, fCentrality, MuPhi, MuCharge, 0.5 };
+ Double_t fillArrayMuRecM[6] = { MuEta, MuPt, fCentrality, MuPhi, MuCharge, MuMSpecies };
+ Double_t fillArrayMuM[6] = { MuMEta, MuMPt, fCentrality, MuMPhi, MuMCharge, MuMSpecies };
+
+ if(fPlotMode==2) { fHMuFM->Fill(fillArrayMuM); fHMuFMrec->Fill(fillArrayMuRecM); }
+ else if(fPlotMode==3) { fHMuPP->Fill(fillArrayMuM); fHMuPPrec->Fill(fillArrayMuRecM); }
+ else if(fPlotMode==4) { fHMuRec->Fill(fillArrayMuRec); }
+ }// end of reconstructed loop
}
-
PostData(1, fOutputList);
return;
}
}
//________________________________________________________________________
-Int_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
+Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
{
- // Applying track cuts for MUON tracks
- Int_t cutnum = 0;
- if(track.ContainTrackerData()) cutnum = 1;
- else return cutnum;
-
- if(track.GetMatchTrigger() > 0) cutnum = 2;
- else return cutnum;
-
- Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
- Double_t eta = track.Eta();
-
- // Eta cut
- if(eta > -4. && eta < -2.5) cutnum = 3;
- else return cutnum;
-
- // Theta cut at absorber end
- if(thetaTrackAbsEnd > 2. && thetaTrackAbsEnd < 10.) cutnum = 4;
- else return cutnum;
-
- // Chi-square cut
- if(track.GetNormalizedChi2() < fChiSquareNormCut) cutnum = 5;
- else return cutnum;
-
- return cutnum;
+ return fMuonTrackCuts->IsSelected(&track);
}
//________________________________________________________________________
-Int_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
+Bool_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
{
- Int_t cutnum = 0;
- if(track.IsMuonTrack()) cutnum = 1;
- else return cutnum;
-
- if (track.GetMatchTrigger() > 0) cutnum = 2;
- else return cutnum;
-
- Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
- Double_t eta = track.Eta();
-
- // Eta cut
- if(eta > -4. && eta < -2.5) cutnum = 3;
- else return cutnum;
-
- // Theta cut at absorber end
- if(thetaTrackAbsEnd > 2. && thetaTrackAbsEnd < 10.) cutnum = 4;
- else return cutnum;
-
- // Chi-square cut
- if(track.Chi2perNDF() < fChiSquareNormCut) cutnum = 5;
- else return cutnum;
-
- return cutnum;
+ return fMuonTrackCuts->IsSelected(&track);
}
//________________________________________________________________________
-void AliMuonEffMC::MDProcess(Bool_t isprimary, Int_t motherpdg, Double_t mcpt, Double_t mcphi, Double_t mceta, Double_t trackpt, Double_t trackphi, Double_t tracketa, Double_t motherpt, Double_t motherphi, Double_t mothereta, Double_t dcavalue)
+Double_t AliMuonEffMC::GetMuonTrackType(AliMCParticle &track)
{
- Int_t motherbin = GetMotherBin(motherpdg);
- Int_t genptbin = -1;
- Int_t recptbin = -1;
-
- if((0. <= mcpt) && (mcpt < 0.5)) genptbin = 0;
- else if((0.5 <= mcpt) && (mcpt < 2.0)) genptbin = 1;
- else genptbin = 2;
+ if(track.GetMother() < fStack->GetNprimary() && track.PdgCode() == 13) return 0.5;
+ else if(track.GetMother() >= fStack->GetNprimary() && track.PdgCode() == 13) return 1.5;
+ else return 2.5;
+}
- if((0. <= trackpt) && (trackpt < 0.5)) recptbin = 0;
- else if((0.5 <= trackpt) && (trackpt < 2.0)) recptbin = 1;
- else recptbin = 2;
- if(isprimary)
- {
- fHMuMotherGenPt[motherbin]->Fill(mcpt, motherpt);
- fHMuMotherRecPt[motherbin]->Fill(trackpt, motherpt);
- fHMuMotherGenPhi[motherbin]->Fill(mcphi, motherphi);
- fHMuMotherRecPhi[motherbin]->Fill(trackphi, motherphi);
- fHMuMotherGenEta[motherbin]->Fill(mceta, mothereta);
- fHMuMotherRecEta[motherbin]->Fill(tracketa, mothereta);
-
- // generated
- fHMuMohterPhiDifGen[motherbin][genptbin]->Fill(deltaphi(motherphi-mcphi));
- fHMuMohterEtaDifGen[motherbin][genptbin]->Fill(mothereta-mceta);
-
- // reconstructed
- fHMuMohterPhiDifRec[motherbin][recptbin]->Fill(deltaphi(motherphi-trackphi));
- fHMuMohterEtaDifRec[motherbin][recptbin]->Fill(mothereta-tracketa);
-
- // DCA
- if(fESD) fHMuDCA[motherbin]->Fill(dcavalue);
- }
- else
- {
- fHMuMotherGenPtSec[motherbin]->Fill(mcpt, motherpt);
- fHMuMotherRecPtSec[motherbin]->Fill(trackpt, motherpt);
- fHMuMotherGenPhiSec[motherbin]->Fill(mcphi, motherphi);
- fHMuMotherRecPhiSec[motherbin]->Fill(trackphi, motherphi);
- fHMuMotherGenEtaSec[motherbin]->Fill(mceta, mothereta);
- fHMuMotherRecEtaSec[motherbin]->Fill(tracketa, mothereta);
-
- // generated
- fHMuMohterPhiDifGenSec[motherbin][genptbin]->Fill(deltaphi(motherphi-mcphi));
- fHMuMohterEtaDifGenSec[motherbin][genptbin]->Fill(mothereta-mceta);
-
- // reconstructed
- fHMuMohterPhiDifRecSec[motherbin][recptbin]->Fill(deltaphi(motherphi-trackphi));
- fHMuMohterEtaDifRecSec[motherbin][recptbin]->Fill(mothereta-tracketa);
-
- // DCA
- if(fESD) fHMuDCASec[motherbin]->Fill(dcavalue);
- }
+//________________________________________________________________________
+Double_t AliMuonEffMC::GetSpecies(Int_t PdgCode)
+{
+ Int_t code = TMath::Abs(PdgCode);
+ if(code==13) return 0.5;
+ else if(code==211) return 1.5;
+ else if(code==321 || code==311) return 2.5;
+ else if(code==411 || code==421) return 3.5;
+ else if(code==511 || code==521) return 4.5;
+ else if(code==213) return 5.5;
+ else if(code==313 || code==323) return 6.5;
+ else if(code==413 || code==423 || code==431) return 7.5;
+ else if(code==513 || code==523 || code==531 || code==533 || code==541 || code==543) return 8.5;
+ else if(code==111) return 9.5;
+ else if(code==113) return 10.5;
+ else if(code==221 || code==331) return 11.5;
+ else if(code==223) return 12.5;
+ else if(code==2212) return 13.5;
+ else if(code==2112) return 14.5;
+ else if(code==1114 || code==2114 || code==2214 || code==2224) return 15.5;
+ else if(code==3112 || code==3114 || code==3212) return 16.5;
+ else if(code>100 && code<1000) return 17.5;
+ else if(code>1000 && code<10000) return 18.5;
+ else return 19.5;
}
//________________________________________________________________________
}
//________________________________________________________________________
-Int_t AliMuonEffMC::GetMotherBin(Int_t motherpdg)
-{
- Int_t motherbin = -1;
-
- if(motherpdg==211) motherbin = 0;
- else if(motherpdg==321) motherbin = 1;
- else if(motherpdg==411 || motherpdg==413 || motherpdg==421 || motherpdg==423 || motherpdg==431 || motherpdg==433 || motherpdg==10413 || motherpdg==10411 || motherpdg==10423 || motherpdg==10421 || motherpdg==10433 || motherpdg==10431 || motherpdg==20413 || motherpdg==415 || motherpdg==20423 || motherpdg==425 || motherpdg==20433 || motherpdg==435) motherbin = 2;
- else motherbin = 3;
-
- return motherbin;
-}
-
-//________________________________________________________________________
-void AliMuonEffMC::FeynmanX()
-{
- for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
- {
- TParticle *particle = fStack->Particle(ipart);
- if(particle->GetFirstMother()==0 || particle->GetFirstMother()==1)
- {
- Int_t pdgcode = particle->GetPdgCode();
- if(pdgcode==2) fHFXu->Fill(TMath::Abs(particle->Pz())/1380.0);
- else if(pdgcode==-2) fHFXantiu->Fill(TMath::Abs(particle->Pz())/1380.0);
- else if(pdgcode==1) fHFXd->Fill(TMath::Abs(particle->Pz())/1380.0);
- else if(pdgcode==-1) fHFXantid->Fill(TMath::Abs(particle->Pz())/1380.0);
- else if(TMath::Abs(pdgcode==21)) fHFXg->Fill(TMath::Abs(particle->Pz())/1380.0);
- else fHFXetc->Fill(TMath::Abs(particle->Pz())/1380.0);
- }
- }
-
- if(fESD)
- {
- for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfMuonTracks(); iTrack++)
- {
- AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
- if(muonTrack)
- {
- if(!IsGoodMUONtrack(*muonTrack)) continue;
- Int_t label = muonTrack->GetLabel();
- if(TMath::Abs(label) > fMC->GetNumberOfTracks()) continue;
-
- TParticle *particle = fStack->Particle(TMath::Abs(label));
- Int_t motherlabel = particle->GetFirstMother();
- if(!(motherlabel==-1) && !(motherlabel==0) && !(motherlabel==1))
- {
- while(1)
- {
- particle = fStack->Particle(motherlabel);
- motherlabel = particle->GetFirstMother();
- if(motherlabel==-1 || motherlabel==0 || motherlabel==1) break;
- }
- }
- if(motherlabel==-1) continue;
-
- Int_t ptbin = -1;
- Int_t etabin = -1;
-
- if(muonTrack->Pt() < 2.0) ptbin=0;
- else ptbin=1;
-
- if(muonTrack->Eta() <= -3.5) etabin=0;
- else if(-3.5<muonTrack->Eta() && muonTrack->Eta()<=-3.0) etabin=1;
- else etabin=2;
-
- if(motherlabel==0) fHFXmuonP[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
- else if(motherlabel==1) fHFXmuonM[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
- }
- }
- }
- else if(fAOD)
- {
- for (Int_t iTrack = 0; iTrack<fAOD->GetNTracks(); iTrack++)
- {
- AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
- if(muonTrack)
- {
- if(!(IsGoodMUONtrack(*muonTrack)) || !(muonTrack->IsMuonTrack())) continue;
- Int_t label = muonTrack->GetLabel();
- if(TMath::Abs(label) > fMC->GetNumberOfTracks()) continue;
-
- TParticle *particle = fStack->Particle(TMath::Abs(label));
- Int_t motherlabel = particle->GetFirstMother();
- if(!(motherlabel==-1) && !(motherlabel==0) && !(motherlabel==1))
- {
- while(1)
- {
- particle = fStack->Particle(motherlabel);
- motherlabel = particle->GetFirstMother();
- if(motherlabel==-1 || motherlabel==0 || motherlabel==1) break;
- }
- }
- if(motherlabel==-1) continue;
-
- Int_t ptbin = -1;
- Int_t etabin = -1;
- if(muonTrack->Pt() < 2.0) ptbin=0;
- else ptbin=1;
-
- if(muonTrack->Eta() <= -3.5) etabin=0;
- else if(-3.5<muonTrack->Eta() && muonTrack->Eta()<=-3.0) etabin=1;
- else etabin=2;
-
- if(motherlabel==0) fHFXmuonP[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
- else if(motherlabel==1) fHFXmuonM[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
- }
- }
- }
-}
-//________________________________________________________________________
-void AliMuonEffMC::ScatFX()
+Int_t AliMuonEffMC::GetFirstMother(Int_t muonlabel)
{
- TParticle *parta = fStack->Particle(6);
- TParticle *partb = fStack->Particle(7);
- if(parta->GetFirstMother()!=-1 || parta->GetFirstMother()!=-1) return;
-
- Double_t xa = TMath::Abs(0.5*(parta->Energy() + partb->Energy() + parta->Pz() + partb->Pz())/1380.0);
- Double_t xb = TMath::Abs(0.5*(parta->Energy() + partb->Energy() - parta->Pz() - partb->Pz())/1380.0);
-
- fHaFx->Fill(xa);
- fHbFx->Fill(xb);
- fHabFxRatio->Fill(xa/xb);
- fHabDeltaFx->Fill(xa-xb);
- fHabRelDeltaFx->Fill((xa-xb)/(xa+xb));
-
- if(fESD)
+ if(fAOD) return -1;
+ else if(fESD)
{
- for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfMuonTracks(); iTrack++)
+ AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(muonlabel);
+ Int_t motherlabel = McParticle->GetMother();
+ Int_t nextmother = 0;
+ if(motherlabel==-1) return -1;
+ else
{
- Int_t MuonPtBin = 0;
- AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
- if(muonTrack)
+ while(motherlabel>-1)
{
- if(!IsGoodMUONtrack(*muonTrack)) continue;
-
- if(muonTrack->Pt() < 1) MuonPtBin = 0;
- else if(muonTrack->Pt() < 2) MuonPtBin = 1;
- else MuonPtBin = 2;
-
- fHabFxMu[MuonPtBin]->Fill(xa, xb);
- fHabFxRatioMu[MuonPtBin]->Fill(xa/xb);
- fHabDeltaFxMu[MuonPtBin]->Fill(xa-xb);
- fHabRelDeltaFxMu[MuonPtBin]->Fill((xa-xb)/(xa+xb));
+ AliMCParticle *MotherParticle = (AliMCParticle*)fMC->GetTrack(motherlabel);
+ nextmother = motherlabel;
+ motherlabel = MotherParticle->GetMother();
}
+ return nextmother;
}
}
-}
-//________________________________________________________________________
-Int_t AliMuonEffMC::GetZVertexBin(Double_t zvertex)
-{
- if(TMath::Abs(zvertex) < 10.0) { return 0; }
- else if(TMath::Abs(zvertex) < 90.0) { return 1; }
- else if(TMath::Abs(zvertex) < 503.0) { return 2; }
- else { return 3; }
+ else return -1;
}
//________________________________________________________________________
-Int_t AliMuonEffMC::GetFirstPrimaryMother(Int_t muonlabel)
+Int_t AliMuonEffMC::GetFirstPPMother(Int_t muonlabel)
{
if(fAOD) return 1;
else if(fESD)
{
AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(muonlabel);
- if(McParticle->GetMother()<fStack->GetNprimary()) return McParticle->GetMother();
+ if(fMC->IsPhysicalPrimary(muonlabel)) return muonlabel;
else
{
Int_t motherlabel = McParticle->GetMother();
while(motherlabel > -1)
{
AliMCParticle *MotherParticle = (AliMCParticle*)fMC->GetTrack(motherlabel);
- if(MotherParticle->GetMother()<fStack->GetNprimary()) break;
+ if(fMC->IsPhysicalPrimary(motherlabel)) break;
else motherlabel = MotherParticle->GetMother();
}
- AliMCParticle *FirstSecondaryMotherParticle = (AliMCParticle*)fMC->GetTrack(motherlabel);
- return FirstSecondaryMotherParticle->GetMother();
- }
- }
- else return -1;
-}
-
-//________________________________________________________________________
-Double_t AliMuonEffMC::GetFirstPrimaryVertex(Int_t muonlabel)
-{
- if(fAOD) return 1.0;
- else if(fESD)
- {
- AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(muonlabel);
- if(McParticle->GetMother()<fStack->GetNprimary()) return McParticle->Zv();
- else
- {
- Int_t motherlabel = McParticle->GetMother();
- while(motherlabel > -1)
- {
- AliMCParticle *MotherParticle = (AliMCParticle*)fMC->GetTrack(motherlabel);
- if(MotherParticle->GetMother()<fStack->GetNprimary()) break;
- else motherlabel = MotherParticle->GetMother();
- }
- AliMCParticle *FirtSecondaryMotherParticle = (AliMCParticle*)fMC->GetTrack(motherlabel);
- return FirtSecondaryMotherParticle->Zv();
+ return motherlabel;
}
}
else return -1;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff