#include "TStyle.h"
//#include "TMCProcess.h"
#include "TArrayI.h"
+#include "TArrayD.h"
#include "TPaveStats.h"
#include "TFitResultPtr.h"
#include "TFile.h"
((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz);
}
- TF1* beautyMuWgt = static_cast<TF1*>(GetWeight("beautyMu"));
-
// Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes // UNCOMMENT TO REJECT PILEUP
// if ( isPileupFromSPD ) return; // UNCOMMENT TO REJECT PILEUP
TObjArray selectedTracks(nTracks);
TArrayI trackSources(nTracks);
- TArrayF trackWgt(nTracks);
+ TArrayD trackWgt(nTracks);
trackWgt.Reset(1.);
for (Int_t itrack = 0; itrack < nTracks; itrack++) {
track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack);
// (before ISR, FSR and kt kick) and the final state one.
// The first muon is of course there only for information and should be rejected.
// The Pythia code for initial state particles is 21
- Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 && AliAnalysisMuonUtility::GetStatusCode(track) != 21 );
+ // When running with POWHEG, Pythia puts the hard process input of POWHEG in the stack
+ // with state 21, and then re-add it to stack before applying ISR, FSR and kt kick.
+ // This muon produces the final state muon, and its status code is 11
+ // To avoid all problems, keep only final state muon (status code <10)
+ // FIXME: is the convention valid for other generators as well?
+ Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 && AliAnalysisMuonUtility::GetStatusCode(track) < 10 );
if ( ! isSelected ) continue;
selectedTracks.AddAt(track,itrack);
trackSources[itrack] = GetParticleType(track);
- if ( trackSources[itrack] == kBeautyMu && beautyMuWgt ) {
- AliVParticle* mcTrack = ( kStepReconstructed ) ? MCEvent()->GetTrack(track->GetLabel()) : track;
- trackWgt[itrack] = beautyMuWgt->Eval(mcTrack->Pt());
+ TObject* wgtObj = GetWeight(fSrcKeys->At(trackSources[itrack])->GetName());
+
+ if ( wgtObj ) {
+ AliVParticle* mcTrack = ( istep == kStepReconstructed ) ? MCEvent()->GetTrack(track->GetLabel()) : track;
+ if ( wgtObj->IsA() == TF1::Class() ) trackWgt[itrack] = static_cast<TF1*>(wgtObj)->Eval(mcTrack->Pt());
+ else if ( wgtObj->IsA()->InheritsFrom(TH1::Class()) ) {
+ TH1* wgtHisto = static_cast<TH1*>(wgtObj);
+ trackWgt[itrack] = wgtHisto->GetBinContent(wgtHisto->GetXaxis()->FindBin(mcTrack->Pt()));
+ }
+ AliDebug(3,Form("Apply weights %s: pt %g gen pt %g weight %g",wgtObj->GetName(),track->Pt(),mcTrack->Pt(),trackWgt[itrack]));
// Int_t iancestor = fUtilityMuonAncestor->GetAncestor(track,MCEvent());
// AliVParticle* motherPart = MCEvent()->GetTrack(iancestor);
// trackWgt[itrack] = beautyMuWgt->GetBinContent(beautyMuWgt->GetXaxis()->FindBin(motherPart->Pt()));
Bool_t isMC = furtherOpt.Contains("MC");
TAxis* srcAxis = gridSparseArray[0]->GetAxis(kHvarMotherType);
- Int_t unIdBin = srcAxis->GetNbins();
- for ( Int_t ibin=1; ibin<=srcAxis->GetNbins(); ibin++ ) {
- currName = srcAxis->GetBinLabel(ibin);
- if ( currName.Contains("Unidentified") ) unIdBin = ibin;
- }
+ Int_t unIdBin = srcAxis->FindBin(fSrcKeys->At(kUnidentified)->GetName());
+ if ( unIdBin < 1 ) unIdBin = srcAxis->GetNbins();
Int_t firstSrcBin = ( isMC ) ? 1 : unIdBin;
Int_t lastSrcBin = ( isMC ) ? srcAxis->GetNbins() - 1 : unIdBin;
}
TH1* histo = gridSparseArray[igrid]->Project(kHvarPt, kHvarEta);
histo->SetName(Form("hPtEta_%s_%s_%s", gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
+ histo->SetDirectory(0);
if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
for ( Int_t iproj=0; iproj<4; ++iproj ) {
histo = gridSparseArray[igrid]->Project(iproj);
histo->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
+ histo->SetDirectory(0);
if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
} // loop on projections
} // loop on grid sparse
axisTitle.ReplaceAll("MuMinus","#mu^{-}");
histo->GetYaxis()->SetTitle(axisTitle.Data());
histo->SetStats(kFALSE);
+ histo->SetDirectory(0);
histoListRatio.Add(histo);
}
evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
-
TString outFilename = Form("/tmp/out%s.root", GetName());
+ printf("\nCreating file %s\n", outFilename.Data());
TFile* outFile = new TFile(outFilename.Data(),"RECREATE");
for ( Int_t icharge=0; icharge<3; icharge++ ) {
histoList[icharge].Write();
}
histoListRatio.Write();
outFile->Close();
- printf("\nCreating file %s\n", outFilename.Data(
- ));
///////////////////
// Vertex method //
sumMothers[kRecoHF].Set(0);
sumMothers[kRecoBkg].Set(0);
sumMothers[kInputHF].Set(3);
- sumMothers[kInputHF][0] = kCharmMu;
- sumMothers[kInputHF][1] = kBeautyMu;
- sumMothers[kInputHF][2] = kQuarkoniumMu;
+
+ sumMothers[kInputHF][0] = srcAxis->FindBin(fSrcKeys->At(kCharmMu)->GetName());
+ sumMothers[kInputHF][1] = srcAxis->FindBin(fSrcKeys->At(kBeautyMu)->GetName());
+ sumMothers[kInputHF][2] = srcAxis->FindBin(fSrcKeys->At(kQuarkoniumMu)->GetName());
sumMothers[kInputDecay].Set(1);
- sumMothers[kInputDecay][0] = kDecayMu;
+ sumMothers[kInputDecay][0] = srcAxis->FindBin(fSrcKeys->At(kDecayMu)->GetName());
sumMothers[kRecoAll].Set(srcAxis->GetNbins());
- for ( Int_t isrc=1; isrc<srcAxis->GetNbins(); ++isrc ) {
- sumMothers[kRecoAll][isrc-1] = isrc;
+ for ( Int_t isrc=0; isrc<srcAxis->GetNbins(); ++isrc ) {
+ sumMothers[kRecoAll][isrc] = isrc+1;
}
meanZ = vtxFit->GetParameter(1);
histo->Divide(eventVertex);
Double_t norm = histo->GetBinContent(histo->FindBin(0.));
histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)");
+ histo->SetTitle(Form("%g < p_{T} (GeV/c) < %g",ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin)));
slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ );