/// Analysis task for single muons in the spectrometer.
/// The output is a list of histograms and CF containers.
/// The macro class can run on AODs or ESDs.
-/// In the latter case a flag can be activated to produce a tree as output.
/// If Monte Carlo information is present, some basics checks are performed.
///
/// \author Diego Stocco
#include "TF1.h"
#include "TStyle.h"
//#include "TMCProcess.h"
-#include "TMinuit.h"
#include "TArrayI.h"
+#include "TArrayD.h"
#include "TPaveStats.h"
+#include "TFitResultPtr.h"
+#include "TFile.h"
+#include "THashList.h"
// STEER includes
#include "AliAODEvent.h"
#include "AliVAnalysisMuon.h"
#include "AliMergeableCollection.h"
#include "AliCounterCollection.h"
+#include "AliMuonEventCuts.h"
#include "AliMuonTrackCuts.h"
+#include "AliAnalysisMuonUtility.h"
/// \cond CLASSIMP
//________________________________________________________________________
AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() :
AliVAnalysisMuon(),
- fMinNvtxContirbutors(0),
- fThetaAbsKeys(0x0)
+ fThetaAbsKeys(0x0),
+ fCutOnDimu(kFALSE)
{
/// Default ctor.
}
//________________________________________________________________________
AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
AliVAnalysisMuon(name, cuts),
- fMinNvtxContirbutors(1),
- fThetaAbsKeys(0x0)
+ fThetaAbsKeys(0x0),
+ fCutOnDimu(kFALSE)
{
//
/// Constructor.
Int_t nPtBins = 80;
Double_t ptMin = 0., ptMax = 80.;
- TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
+ TString ptName("Pt"), ptTitle("p_{T}"), ptUnits("GeV/c");
Int_t nEtaBins = 25;
Double_t etaMin = -4.5, etaMax = -2.;
TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
Int_t nPhiBins = 36;
- Double_t phiMin = 0.; Double_t phiMax = 2*TMath::Pi();
+ Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
Int_t nChargeBins = 2;
AddObjectToCollection(cfContainer, kTrackContainer);
+ histoName = "hRecoDimu";
+ TH2* histoDimu = new TH2F(histoName.Data(), histoName.Data(), 24, 0., 120., 30, 0., 150.);
+ histoDimu->SetXTitle("min. p_{T}^{#mu} (GeV/c)");
+ histoDimu->SetYTitle("M_{#mu#mu} (GeV/c^{2})");
+ AddObjectToCollection(histoDimu, kNobjectTypes);
+
+ histoName = "hGeneratedZ";
+ 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->SetXTitle("#eta_{#mu-}");
+ histoDimu->SetYTitle("#eta_{#mu+}");
+ AddObjectToCollection(histoDimu, kNobjectTypes+2);
+
fMuonTrackCuts->Print("mask");
+
+ AliInfo(Form("Apply cut on dimuon (60<M_mumu<120 GeV/c^2) to reject Z contribution: %i", fCutOnDimu));
}
//________________________________________________________________________
/// Fill output objects
//
- AliVVertex* primaryVertex = ( fAODEvent ) ? (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() : (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
- if ( primaryVertex->GetNContributors() < fMinNvtxContirbutors ) return;
-
- Double_t ipVz = primaryVertex->GetZ();
- Double_t ipVzMC = 0;
- if ( IsMC() ) {
- if ( fMCEvent ) ipVzMC = fMCEvent->GetPrimaryVertex()->GetZ();
- else if ( fAODEvent ) {
- AliAODMCHeader* aodMCHeader = (AliAODMCHeader *)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName());
- if ( aodMCHeader ) ipVzMC = aodMCHeader->GetVtxZ();
- }
- }
+ Double_t ipVz = AliAnalysisMuonUtility::GetVertexSPD(InputEvent())->GetZ();
+ Double_t ipVzMC = MCEvent() ? AliAnalysisMuonUtility::GetMCVertexZ(InputEvent(),MCEvent()) : 0.;
for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
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
- if ( isPileupFromSPD ) return;
+
+ // 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
Double_t containerInput[kNvars];
AliVParticle* track = 0x0;
- for ( Int_t istep = 0; istep<2; ++istep ) {
- Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
+ Int_t nSteps = MCEvent() ? 2 : 1;
+ for ( Int_t istep = 0; istep<nSteps; ++istep ) {
+ Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : MCEvent()->GetNumberOfTracks();
+
+ TObjArray selectedTracks(nTracks);
+ TArrayI trackSources(nTracks);
+ TArrayD trackWgt(nTracks);
+ trackWgt.Reset(1.);
for (Int_t itrack = 0; itrack < nTracks; itrack++) {
- track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
+ track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack);
- Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
+ // In case of MC we usually ask that the particle is a muon
+ // However, in W or Z simulations, Pythia stores both the initial muon
+ // (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 );
if ( ! isSelected ) continue;
- // In W simulations with Pythia, sometimes muon is stored twice.
- // Remove muon in case it has another muon as daugther
- if ( istep == kStepGeneratedMC ) {
- Int_t firstDaughter = GetDaughterIndex(track, 0);
- if ( firstDaughter >= 0 ) {
- Bool_t hasMuonDaughter = kFALSE;
- Int_t lastDaughter = GetDaughterIndex(track, 1);
- for ( Int_t idaugh=firstDaughter; idaugh<=lastDaughter; idaugh++ ) {
- AliVParticle* currTrack = GetMCTrack(idaugh);
- if ( currTrack->PdgCode() == track->PdgCode() ) {
- hasMuonDaughter = kTRUE;
- break;
- }
- }
- if ( hasMuonDaughter ) {
- AliDebug(1, Form("Current muon (%i) has muon daughter: rejecting it", itrack));
- 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
+ TArrayI rejectTrack(nTracks);
+ for ( Int_t itrack=0; itrack<nTracks; itrack++) {
+ track = static_cast<AliVParticle*>(selectedTracks.At(itrack));
+ if ( ! track ) continue;
- Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
+ // 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 ( 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);
+
+ 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;
if ( istep == kStepReconstructed ) {
- Double_t rAbsEnd = ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
- thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
+ thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
}
else {
thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg();
containerInput[kHvarVz] = ( istep == kStepReconstructed ) ? ipVz : ipVzMC;
containerInput[kHvarCharge] = track->Charge()/3.;
containerInput[kHvarThetaAbs] = (Double_t)thetaAbsBin;
- containerInput[kHvarMotherType] = (Double_t)trackSrc;
+ containerInput[kHvarMotherType] = (Double_t)trackSources[itrack];
for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
- if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
- ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep);
+ if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
+ ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep,trackWgt[itrack]);
} // loop on selected trigger classes
} // loop on tracks
} // loop on container steps
//
/// Draw some histograms at the end.
//
-
+
AliVAnalysisMuon::Terminate("");
-
+
if ( ! fMergeableCollection ) return;
TString physSel = fTerminateOptions->At(0)->GetName();
TString trigClassName = fTerminateOptions->At(1)->GetName();
TString centralityRange = fTerminateOptions->At(2)->GetName();
TString furtherOpt = fTerminateOptions->At(3)->GetName();
+
+ TString minBiasTrig = "";
+ TObjArray* optArr = furtherOpt.Tokenize(" ");
+ TString currName = "";
+ for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
+ currName = optArr->At(iopt)->GetName();
+ if ( currName.Contains("-B-") || currName.Contains("ANY") ) minBiasTrig = currName;
+ }
+ delete optArr;
+
furtherOpt.ToUpper();
+ Bool_t plotChargeAsymmetry = furtherOpt.Contains("ASYM");
+
+ AliCFContainer* inContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") );
+ if ( ! inContainer ) return;
+
+ AliCFContainer* cfContainer = inContainer;
- AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") );
- if ( ! cfContainer ) return;
+ if ( ! furtherOpt.Contains("GENINTRIGCLASS") && trigClassName != "ANY" ) {
+ // The output container contains both the reconstructed and (in MC)
+ // the generated muons in a specific trigger class.
+ // Since the trigger pt level of the track is matched to the trigger class,
+ // analyzing the MUHigh trigger (for example) is equivalent of analyzing
+ // Hpt tracks.
+ // However, in this way, the generated muons distribution depend
+ // on a condition on the reconstructed track.
+ // If we calculate the Acc.xEff. as the ratio of reconstructed and
+ // generated muons IN A TRIGGER CLASS, we are biasing the final result.
+ // The correct value of Acc.xEff. is hence obtained as the distribution
+ // of reconstructed tracks in a muon trigger class, divided by the
+ // total number of generated class (which is in the "class" ANY).
+ // The following code sets the generated muons as the one in the class ANY.
+ // The feature is the default. If you want the generated muons in the same
+ // trigger class as the generated tracks, use option "GENINTRIGCLASS"
+ AliCFContainer* fullMCcontainer = static_cast<AliCFContainer*> ( GetSum(physSel,"ANY",centralityRange,"SingleMuContainer") );
+ if ( fullMCcontainer ) {
+ cfContainer = static_cast<AliCFContainer*>(cfContainer->Clone("SingleMuContainerCombo"));
+ cfContainer->SetGrid(kStepGeneratedMC,fullMCcontainer->GetGrid(kStepGeneratedMC));
+ }
+ }
AliCFEffGrid* effSparse = new AliCFEffGrid(Form("eff%s", cfContainer->GetName()),Form("Efficiency %s", cfContainer->GetTitle()),*cfContainer);
effSparse->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC);
AliCFGridSparse* gridSparseArray[3] = {effSparse->GetNum(), effSparse->GetDen(), effSparse};
TString gridSparseName[3] = {cfContainer->GetStepTitle(kStepReconstructed), cfContainer->GetStepTitle(kStepGeneratedMC), "Efficiency"};
-
- Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
-// TString allSrcNames = "";
-// for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
-// if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" ");
-// allSrcNames += fSrcKeys->At(isrc)->GetName();
-// }
-
+
+ Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
+ // TString allSrcNames = "";
+ // for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
+ // if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" ");
+ // allSrcNames += fSrcKeys->At(isrc)->GetName();
+ // }
+
TCanvas* can = 0x0;
Int_t xshift = 100;
Int_t yshift = 100;
Int_t igroup2 = 0;
Bool_t isMC = furtherOpt.Contains("MC");
- Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
- Int_t lastSrc = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
- if ( ! isMC ) srcColors[kUnidentified] = 1;
-
- TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
+
+ TAxis* srcAxis = gridSparseArray[0]->GetAxis(kHvarMotherType);
+ 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;
+ if ( ! isMC ) srcColors[unIdBin-1] = 1;
+
////////////////
// Kinematics //
////////////////
- TCanvas* canKine[3] = {0x0, 0x0, 0x0};
- TLegend* legKine[3] = {0x0, 0x0, 0x0};
- for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
- for ( Int_t icharge=0; icharge<2; ++icharge ) {
+ TString chargeNames[3] = {fChargeKeys->At(0)->GetName(), fChargeKeys->At(1)->GetName(), "Total"};
+ THashList histoList[3];
+ for ( Int_t icharge=0; icharge<3; icharge++ ) {
+ histoList[icharge].SetName(chargeNames[icharge].Data());
+ }
+ for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) {
+ for ( Int_t icharge=0; icharge<3; ++icharge ) {
+ Int_t icharge1 = ( icharge < 2 ) ? icharge : 0;
+ Int_t icharge2 = ( icharge < 2 ) ? icharge : 1;
for ( Int_t igrid=0; igrid<3; ++igrid ) {
if ( gridSparseArray[igrid]->GetEntries() == 0. ) break;
if ( gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) {
SetSparseRange(gridSparseArray[igrid], kHvarEta, "", -3.999, -2.501);
- SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
- SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge+1, icharge+1, "USEBIN");
- }
- if ( ! canKine[igrid] ) {
- igroup1++;
- igroup2 = 0;
- currName = Form("%s_proj_%s", GetName(), gridSparseName[igrid].Data());
- canKine[igrid] = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
- canKine[igrid]->Divide(2,2);
- legKine[igrid] = new TLegend(0.6, 0.6, 0.8, 0.8);
- igroup2++;
+ SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc, isrc, "USEBIN");
+ SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge1+1, icharge2+1, "USEBIN");
}
+ TH1* histo = gridSparseArray[igrid]->Project(kHvarPt, kHvarEta);
+ histo->SetName(Form("hPtEta_%s_%s_%s", gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
+ if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
for ( Int_t iproj=0; iproj<4; ++iproj ) {
- canKine[igrid]->cd(iproj+1);
- if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy();
- TH1* projHisto = gridSparseArray[igrid]->Project(iproj);
- projHisto->SetName(Form("proj%i_%s_src%i_charge%i", iproj, gridSparseName[igrid].Data(), isrc, icharge));
- if ( projHisto->GetEntries() == 0 ) continue;
- Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 );
- drawOpt = isFirst ? "e" : "esames";
- //if ( isrc == kUnidentified && ! drawOpt.Contains("same") ) isMC = kFALSE;
- //if ( ! isMC ) srcColors[kUnidentified] = 1;
- projHisto->SetLineColor(srcColors[isrc]);
- projHisto->SetMarkerColor(srcColors[isrc]);
- projHisto->SetMarkerStyle(20+4*icharge);
- projHisto->Draw(drawOpt.Data());
- gPad->Update();
- TPaveStats* paveStats = (TPaveStats*)projHisto->FindObject("stats");
- if ( paveStats ) paveStats->SetTextColor(srcColors[isrc]);
- if ( iproj == 0 ) {
- TString legEntry = fChargeKeys->At(icharge)->GetName();
- if ( isMC ) legEntry += Form(" %s", fSrcKeys->At(isrc)->GetName());
- legKine[igrid]->AddEntry(projHisto,legEntry.Data(), "lp");
- }
- } // loop on grid sparse
- } // loop on projections
- } // loop on mu charge
+ histo = gridSparseArray[igrid]->Project(iproj);
+ histo->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
+ if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
+ } // loop on projections
+ } // loop on grid sparse
+ } // loop on charge
} // loop on track sources
+ // Get charge asymmetry or mu+/mu-
+ THashList histoListRatio;
+ TString basePlotName = plotChargeAsymmetry ? "ChargeAsymmetry" : "ChargeRatio";
+ histoListRatio.SetName(basePlotName.Data());
+ Int_t baseCharge = 1;
+ Int_t auxCharge = 1-baseCharge;
+ for ( Int_t ihisto=0; ihisto<histoList[baseCharge].GetEntries(); ihisto++ ) {
+ TObject* obj = histoList[baseCharge].At(ihisto);
+ TString histoName = obj->GetName();
+ if ( histoName.Contains(gridSparseName[2].Data()) ) continue;
+ TString searchName = histoName;
+ searchName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(), fChargeKeys->At(auxCharge)->GetName());
+ TH1* auxHisto = static_cast<TH1*> (histoList[auxCharge].FindObject(searchName.Data()));
+ if ( ! auxHisto ) continue;
+ histoName.ReplaceAll(fChargeKeys->At(baseCharge)->GetName(),basePlotName.Data());
+ TH1* histo = static_cast<TH1*> (obj->Clone(histoName.Data()));
+ if ( plotChargeAsymmetry ) {
+ histo->Add(auxHisto, -1.);
+ // h2 + h1 = 2xh2 + (h1-h2)
+ auxHisto->Add(auxHisto, histo, 2.);
+ }
+ histo->Divide(auxHisto);
+ TString axisTitle = plotChargeAsymmetry ? Form("(%s - %s) / (%s + %s)", fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName(), fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName()) : Form("%s / %s", fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName());
+ axisTitle.ReplaceAll("MuPlus","#mu^{+}");
+ axisTitle.ReplaceAll("MuMinus","#mu^{-}");
+ histo->GetYaxis()->SetTitle(axisTitle.Data());
+ histo->SetStats(kFALSE);
+ histoListRatio.Add(histo);
+ }
+
+ // Plot kinematics
+ TString histoName = "", drawOpt = "";
+ for ( Int_t itype=0; itype<3; itype++ ) {
+ THashList* currList = 0x0;
+ Int_t nCharges = 1;
+ if ( itype == 1 ) currList = &histoListRatio;
+ else if ( itype == 2 ) currList = &histoList[2];
+ else nCharges = 2;
+ for ( Int_t igrid=0; igrid<3; ++igrid ) {
+ igroup1 = igrid;
+ TCanvas* canKine = 0x0;
+ TLegend* legKine = 0x0;
+ for ( Int_t iproj=0; iproj<4; ++iproj ) {
+ for ( Int_t isrc = firstSrcBin; isrc <= lastSrcBin; ++isrc ) {
+ for ( Int_t icharge=0; icharge<nCharges; ++icharge ) {
+ if ( itype == 0 ) currList = &histoList[icharge];
+ histoName = Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), currList->GetName());
+ TH1* histo = static_cast<TH1*>(currList->FindObject(histoName.Data()));
+ if ( ! histo ) continue;
+ if ( ! canKine ) {
+ igroup2 = itype;
+ igroup1 = igrid;
+ currName = Form("%s_%s_%s", GetName(), currList->GetName(), gridSparseName[igrid].Data());
+ canKine = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+ canKine->Divide(2,2);
+ legKine = new TLegend(0.6, 0.6, 0.8, 0.8);
+ }
+ canKine->cd(iproj+1);
+ if ( itype != 1 ) {
+ if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy();
+ }
+ Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 );
+ drawOpt = isFirst ? "e" : "esames";
+ histo->SetLineColor(srcColors[isrc-1]);
+ histo->SetMarkerColor(srcColors[isrc-1]);
+ histo->SetMarkerStyle(20+4*icharge);
+ histo->Draw(drawOpt.Data());
+ TPaveStats* paveStats = (TPaveStats*)histo->FindObject("stats");
+ if ( paveStats ) paveStats->SetTextColor(srcColors[isrc-1]);
+ if ( iproj == 0 ) {
+ TString legEntry = ( itype == 0 ) ? fChargeKeys->At(icharge)->GetName() : "";
+ if ( isMC ) legEntry += Form(" %s", srcAxis->GetBinLabel(isrc));
+ if ( ! legEntry.IsNull() ) legKine->AddEntry(histo,legEntry.Data(), "lp");
+ }
+ } // loop on mu charge
+ } // loop on track sources
+ } // loop on projections
+ if ( legKine && legKine->GetNRows() > 0 ) {
+ canKine->cd(1);
+ legKine->Draw("same");
+ }
+ } // loop on grid sparse
+ } // loop on types
+
for ( Int_t igrid=0; igrid<3; igrid++ ) {
- if ( ! canKine[igrid] ) continue;
- canKine[igrid]->cd(1);
- legKine[igrid]->Draw("same");
if ( gridSparseArray[igrid]->IsA() == AliCFEffGrid::Class() ) continue;
SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", 1, gridSparseArray[igrid]->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range
} // loop on container steps
-
//////////////////////
// Event statistics //
//////////////////////
fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
+ TString outFilename = Form("/tmp/out%s.root", GetName());
+ 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 = kUnidentified, lastMother = kUnidentified;
igroup1++;
- TH1* eventVertex = (TH1*)GetSum(physSel, "CINT7-I-NOPF-ALLNOTRD", centralityRange, "hIpVtx");
+ TH1* eventVertex = (TH1*)GetSum(physSel, minBiasTrig, centralityRange, "hIpVtx");
if ( ! eventVertex ) return;
Double_t minZ = -9.99, maxZ = 9.99;
Double_t meanZ = 0., sigmaZ = 4.;
Double_t nSigma = 2.;
- TString fitOpt = "R0";
+ TString fitOpt = "R0S";
Bool_t fixFitRange = kFALSE;
TString fitFormula = Form("[0]+[1]*(x+[2])");
-
- // Get vertex shape
+
+ // Get vertex shape
if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2();
Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow
printf("Event vertex integral %.0f\n\n", eventVtxIntegral);
can->SetLogy();
eventVertex->Draw();
vtxFit->Draw("same");
-
+
enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos};
TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"};
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[kRecoAll].Set(kNtrackSources);
- for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
- sumMothers[kRecoAll][isrc] = isrc;
+ sumMothers[kInputDecay][0] = srcAxis->FindBin(fSrcKeys->At(kDecayMu)->GetName());
+ sumMothers[kRecoAll].Set(srcAxis->GetNbins());
+ for ( Int_t isrc=0; isrc<srcAxis->GetNbins(); ++isrc ) {
+ sumMothers[kRecoAll][isrc] = isrc+1;
}
meanZ = vtxFit->GetParameter(1);
const Double_t kFreePath = 153.; // 150.; // 130.; // cm
//fitFunc->SetParameters(0.,1.);
fitFunc->FixParameter(2, kFreePath);
-
+
AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed);
TAxis* ptAxis = gridSparse->GetAxis(kHvarPt);
Double_t slope = 0.;
Double_t limitNorm = 0., limitSlope = 0.;
- Int_t firstPtBin = 0, lastPtBin = 0;
-
+ Int_t firstPtBin = 0, lastPtBin = 0;
+
+ gStyle->SetOptFit(1111);
+
for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
igroup2++;
SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN");
recoHisto[ireco]->Reset();
recoHisto[ireco]->Sumw2();
for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) {
- SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc]+1, sumMothers[ireco][isrc]+1, "USEBIN");
+ SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc], sumMothers[ireco][isrc], "USEBIN");
TH1* auxHisto = gridSparse->Project(kHvarPt);
recoHisto[ireco]->Add(auxHisto);
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 ) {
firstPtBin = ibinpt;
lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt;
SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN");
TH1* histo = gridSparse->Project(kHvarVz);
histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt));
- if ( histo->GetEntries() < 100. ) break;
+ if ( histo->Integral() < 100. ) break;
printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries());
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)");
- slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
+ 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 );
if ( slope < 0. ) slope = norm/kFreePath;
fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK
printf("Norm 0. < %f < %f slope 0. < %f < %f\n", norm, limitNorm, slope, limitSlope);
}
- histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit);
+ TFitResultPtr fitRes = histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit);
- if ( gMinuit->fCstatu.Contains("CONVERGED") &&
- fitFunc->GetParameter(0) > 0. &&
+ // if ( gMinuit->fCstatu.Contains("CONVERGED") &&
+ if ( ((Int_t)fitRes) == 0 &&
+ fitFunc->GetParameter(0) > 0. &&
fitFunc->GetParameter(1) > 0. )
break;
else if ( furtherOpt.Contains("REFIT") ) printf("Re-fit with limits\n");
fitFunc->DrawCopy("same");
currDraw++;
} // loop on pt bins
- SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
- TH1* totalPtDistrib = gridSparse->Project(kHvarPt);
- totalPtDistrib->SetName(Form("totalPtDistrib_%s", fThetaAbsKeys->GetName()));
+ 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