#include "TStyle.h"
//#include "TMCProcess.h"
#include "TArrayI.h"
+#include "TArrayD.h"
#include "TPaveStats.h"
#include "TFitResultPtr.h"
#include "TFile.h"
AddObjectToCollection(histoDimu, kNobjectTypes);
histoName = "hGeneratedZ";
- AddObjectToCollection(histoDimu->Clone(histoName.Data()), kNobjectTypes+1);
+ TH2* histoGenZ = static_cast<TH2*>(histoDimu->Clone(histoName.Data()));
+ histoGenZ->SetTitle(histoName.Data());
+ AddObjectToCollection(histoGenZ, kNobjectTypes+1);
histoName = "hZmuEtaCorr";
- histoDimu = new TH2F(histoName.Data(), histoName.Data(), 160, -8., 8., 160,-8., 8.);
+ histoDimu = new TH2F(histoName.Data(), histoName.Data(), 160, -8., 8., 160, -8., 8.);
histoDimu->SetXTitle("#eta_{#mu-}");
histoDimu->SetYTitle("#eta_{#mu+}");
AddObjectToCollection(histoDimu, kNobjectTypes+2);
TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz);
}
-
+
// 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);
+ 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);
+ 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()));
+ }
} // loop on tracks
// Loop on selected tracks
track = static_cast<AliVParticle*>(selectedTracks.At(itrack));
if ( ! track ) continue;
- if ( istep == kStepGeneratedMC || fCutOnDimu ) {
- // Check dimuons
- for ( Int_t jtrack=itrack+1; jtrack<nTracks; jtrack++ ) {
- AliVParticle* auxTrack = static_cast<AliVParticle*>(selectedTracks.At(jtrack));
- if ( ! auxTrack ) continue;
- if ( track->Charge() * auxTrack->Charge() >= 0 ) continue;
+ // Check dimuons
+ for ( Int_t jtrack=itrack+1; jtrack<nTracks; jtrack++ ) {
+ AliVParticle* auxTrack = static_cast<AliVParticle*>(selectedTracks.At(jtrack));
+ if ( ! auxTrack ) continue;
+ if ( track->Charge() * auxTrack->Charge() >= 0 ) continue;
- TLorentzVector dimuPair = AliAnalysisMuonUtility::GetTrackPair(track,auxTrack);
- Double_t ptMin = TMath::Min(track->Pt(),auxTrack->Pt());
- Double_t invMass = dimuPair.M();
- if ( invMass > 60. && invMass < 120. ) {
- rejectTrack[itrack] = 1;
- rejectTrack[jtrack] = 1;
- }
- if ( istep == kStepReconstructed ) {
- for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
- TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
- if ( ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) || ! fMuonTrackCuts->TrackPtCutMatchTrigClass(auxTrack, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
- ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hRecoDimu"))->Fill(ptMin,invMass);
- } // loop on triggers
- }
- else {
- if ( trackSources[itrack] == kZbosonMu && trackSources[jtrack] == kZbosonMu ) {
- Bool_t isAccepted = kTRUE;
- AliVParticle* muDaughter[2] = {0x0, 0x0};
- for ( Int_t imu=0; imu<2; imu++ ) {
- AliVParticle* currPart = ( imu == 0 ) ? track : auxTrack;
- if ( currPart->Charge() < 0. ) muDaughter[0] = currPart;
- else muDaughter[1] = currPart;
- if ( currPart->Eta() < -4.5 || currPart->Eta() > -2. ) {
- isAccepted = kFALSE;
- }
- } // loop on muons in the pair
+ TLorentzVector dimuPair = AliAnalysisMuonUtility::GetTrackPair(track,auxTrack);
+ Double_t ptMin = TMath::Min(track->Pt(),auxTrack->Pt());
+ Double_t invMass = dimuPair.M();
+ if ( fCutOnDimu && invMass > 60. && invMass < 120. ) {
+ rejectTrack[itrack] = 1;
+ rejectTrack[jtrack] = 1;
+ }
+ Double_t dimuWgt = trackWgt[itrack] * trackWgt[jtrack];
+ if ( istep == kStepReconstructed ) {
+ for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
+ TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
+ if ( ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) || ! fMuonTrackCuts->TrackPtCutMatchTrigClass(auxTrack, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
+ ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hRecoDimu"))->Fill(ptMin,invMass,dimuWgt);
+ } // loop on triggers
+ }
+ else {
+ if ( trackSources[itrack] == kZbosonMu && trackSources[jtrack] == kZbosonMu ) {
+ Bool_t isAccepted = kTRUE;
+ AliVParticle* muDaughter[2] = {0x0, 0x0};
+ for ( Int_t imu=0; imu<2; imu++ ) {
+ AliVParticle* currPart = ( imu == 0 ) ? track : auxTrack;
+ if ( currPart->Charge() < 0. ) muDaughter[0] = currPart;
+ else muDaughter[1] = currPart;
+ if ( currPart->Eta() < -4.5 || currPart->Eta() > -2. ) {
+ isAccepted = kFALSE;
+ }
+ } // loop on muons in the pair
- Double_t pairRapidity = dimuPair.Rapidity();
- if ( pairRapidity < -4. || pairRapidity > -2.5 ) isAccepted = kFALSE;
- //printf("Rapidity Z %g pair %g\n",track->Y(), pairRapidity);
+ Double_t pairRapidity = dimuPair.Rapidity();
+ if ( pairRapidity < -4. || pairRapidity > -2.5 ) isAccepted = kFALSE;
+ //printf("Rapidity Z %g pair %g\n",track->Y(), pairRapidity);
- for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
- TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
- if ( isAccepted ) ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hGeneratedZ"))->Fill(ptMin,invMass);
- ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hZmuEtaCorr"))->Fill(muDaughter[0]->Eta(),muDaughter[1]->Eta());
- } // loop on selected trig
- } // both muons from Z
- } // kStepGeneratedMC
- } // loop on auxiliary tracks
- } // apply cut on dimu
+ for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
+ TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
+ if ( isAccepted ) ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hGeneratedZ"))->Fill(ptMin,invMass,dimuWgt);
+ ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, "hZmuEtaCorr"))->Fill(muDaughter[0]->Eta(),muDaughter[1]->Eta(),dimuWgt);
+ } // loop on selected trig
+ } // both muons from Z
+ } // kStepGeneratedMC
+ } // loop on auxiliary tracks
if ( rejectTrack[itrack] > 0 ) continue;
Double_t thetaAbsEndDeg = 0;
for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
- ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep);
+ ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep,trackWgt[itrack]);
} // loop on selected trigger classes
} // loop on tracks
} // loop on container steps
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 //
///////////////////
if ( ! furtherOpt.Contains("VERTEX") ) return;
- Int_t firstMother = ( isMC ) ? 0 : kUnidentified;
- Int_t lastMother = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
igroup1++;
TH1* eventVertex = (TH1*)GetSum(physSel, minBiasTrig, centralityRange, "hIpVtx");
if ( ! eventVertex ) return;
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);
delete auxHisto;
}
}
- SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
+ SetSparseRange(gridSparse, kHvarMotherType, "", firstPtBin, lastSrcBin, "USEBIN");
Int_t currDraw = 0;
for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) {
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 );
fitFunc->DrawCopy("same");
currDraw++;
} // loop on pt bins
- SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
+ SetSparseRange(gridSparse, kHvarMotherType, "", firstSrcBin, lastSrcBin, "USEBIN");
currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName());
can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600);
TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff