]>
Commit | Line | Data |
---|---|---|
aad6618e | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
27de2dfb | 16 | /* $Id$ */ |
17 | ||
662e37fe | 18 | //----------------------------------------------------------------------------- |
19 | /// \class AliAnalysisTaskSingleMu | |
20 | /// Analysis task for single muons in the spectrometer. | |
a07db2fb | 21 | /// The output is a list of histograms and CF containers. |
b201705a | 22 | /// The macro class can run on AODs or ESDs. |
9728bcfd | 23 | /// If Monte Carlo information is present, some basics checks are performed. |
662e37fe | 24 | /// |
25 | /// \author Diego Stocco | |
26 | //----------------------------------------------------------------------------- | |
27 | ||
aad6618e | 28 | #define AliAnalysisTaskSingleMu_cxx |
29 | ||
a07db2fb | 30 | #include "AliAnalysisTaskSingleMu.h" |
31 | ||
aad6618e | 32 | // ROOT includes |
aad6618e | 33 | #include "TROOT.h" |
9728bcfd | 34 | #include "TH1.h" |
35 | #include "TH2.h" | |
9728bcfd | 36 | #include "TAxis.h" |
662e37fe | 37 | #include "TCanvas.h" |
a07db2fb | 38 | #include "TLegend.h" |
9728bcfd | 39 | #include "TMath.h" |
589e3f71 | 40 | #include "TObjString.h" |
55065f3f | 41 | #include "TObjArray.h" |
a07db2fb | 42 | #include "TF1.h" |
43 | #include "TStyle.h" | |
44 | //#include "TMCProcess.h" | |
a07db2fb | 45 | #include "TArrayI.h" |
c6e0141f | 46 | #include "TPaveStats.h" |
12e33589 | 47 | #include "TFitResultPtr.h" |
aad6618e | 48 | |
49 | // STEER includes | |
aad6618e | 50 | #include "AliAODEvent.h" |
51 | #include "AliAODTrack.h" | |
a07db2fb | 52 | #include "AliAODMCParticle.h" |
75008b31 | 53 | #include "AliMCEvent.h" |
9728bcfd | 54 | #include "AliMCParticle.h" |
b201705a | 55 | #include "AliESDEvent.h" |
56 | #include "AliESDMuonTrack.h" | |
a07db2fb | 57 | #include "AliVHeader.h" |
58 | #include "AliAODMCHeader.h" | |
59 | #include "AliStack.h" | |
9728bcfd | 60 | |
589e3f71 | 61 | // ANALYSIS includes |
62 | #include "AliAnalysisManager.h" | |
589e3f71 | 63 | |
64 | // CORRFW includes | |
589e3f71 | 65 | #include "AliCFContainer.h" |
66 | #include "AliCFGridSparse.h" | |
c6e0141f | 67 | #include "AliCFEffGrid.h" |
589e3f71 | 68 | |
c6e0141f | 69 | // PWG includes |
a07db2fb | 70 | #include "AliVAnalysisMuon.h" |
71 | #include "AliMergeableCollection.h" | |
72 | #include "AliCounterCollection.h" | |
c6e0141f | 73 | #include "AliMuonTrackCuts.h" |
a07db2fb | 74 | |
aad6618e | 75 | |
662e37fe | 76 | /// \cond CLASSIMP |
77 | ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context | |
78 | /// \endcond | |
aad6618e | 79 | |
9728bcfd | 80 | |
aad6618e | 81 | //________________________________________________________________________ |
a07db2fb | 82 | AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() : |
83 | AliVAnalysisMuon(), | |
84 | fThetaAbsKeys(0x0) | |
aad6618e | 85 | { |
a07db2fb | 86 | /// Default ctor. |
aad6618e | 87 | } |
88 | ||
9728bcfd | 89 | //________________________________________________________________________ |
a07db2fb | 90 | AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) : |
91 | AliVAnalysisMuon(name, cuts), | |
92 | fThetaAbsKeys(0x0) | |
aad6618e | 93 | { |
94 | // | |
a07db2fb | 95 | /// Constructor. |
9728bcfd | 96 | // |
a07db2fb | 97 | TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310"; |
98 | fThetaAbsKeys = thetaAbsKeys.Tokenize(" "); | |
589e3f71 | 99 | } |
100 | ||
101 | ||
a07db2fb | 102 | //________________________________________________________________________ |
103 | AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu() | |
589e3f71 | 104 | { |
105 | // | |
a07db2fb | 106 | /// Destructor |
589e3f71 | 107 | // |
108 | ||
a07db2fb | 109 | delete fThetaAbsKeys; |
aad6618e | 110 | } |
111 | ||
112 | //___________________________________________________________________________ | |
a07db2fb | 113 | void AliAnalysisTaskSingleMu::MyUserCreateOutputObjects() |
aad6618e | 114 | { |
589e3f71 | 115 | |
a07db2fb | 116 | TH1* histo = 0x0; |
117 | TString histoName = "", histoTitle = ""; | |
9728bcfd | 118 | |
589e3f71 | 119 | Int_t nVzBins = 40; |
120 | Double_t vzMin = -20., vzMax = 20.; | |
a07db2fb | 121 | TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm"); |
9728bcfd | 122 | |
a07db2fb | 123 | histoName = "hIpVtx"; |
124 | histo = new TH1F(histoName.Data(), histoName.Data(), nVzBins, vzMin, vzMax); | |
125 | histo->SetXTitle("v_{z} (cm)"); | |
126 | AddObjectToCollection(histo, kIPVz); | |
589e3f71 | 127 | |
a07db2fb | 128 | |
129 | Int_t nPtBins = 80; | |
130 | Double_t ptMin = 0., ptMax = 80.; | |
131 | TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c"); | |
132 | ||
133 | Int_t nEtaBins = 25; | |
134 | Double_t etaMin = -4.5, etaMax = -2.; | |
135 | TString etaName("Eta"), etaTitle("#eta"), etaUnits(""); | |
136 | ||
137 | Int_t nPhiBins = 36; | |
12e33589 | 138 | Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi(); |
a07db2fb | 139 | TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad"); |
140 | ||
589e3f71 | 141 | Int_t nChargeBins = 2; |
142 | Double_t chargeMin = -2, chargeMax = 2.; | |
143 | TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e"); | |
9728bcfd | 144 | |
a07db2fb | 145 | Int_t nThetaAbsEndBins = 2; |
146 | Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 1.5; | |
147 | TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u."); | |
148 | ||
589e3f71 | 149 | Int_t nMotherTypeBins = kNtrackSources; |
150 | Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5; | |
151 | TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits(""); | |
a07db2fb | 152 | |
153 | Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nVzBins, nChargeBins, nThetaAbsEndBins, nMotherTypeBins}; | |
154 | Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, vzMin, chargeMin, thetaAbsEndMin, motherTypeMin}; | |
155 | Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, vzMax, chargeMax, thetaAbsEndMax, motherTypeMax}; | |
156 | TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, vzTitle, chargeTitle, thetaAbsEndTitle, motherTypeTitle}; | |
157 | TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, vzUnits, chargeUnits, thetaAbsEndUnits, motherTypeUnits}; | |
9728bcfd | 158 | |
a07db2fb | 159 | AliCFContainer* cfContainer = new AliCFContainer("SingleMuContainer","Container for tracks",kNsteps,kNvars,nbins); |
4ab8d5a6 | 160 | |
589e3f71 | 161 | for ( Int_t idim = 0; idim<kNvars; idim++){ |
162 | histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data()); | |
163 | histoTitle.ReplaceAll("()",""); | |
a07db2fb | 164 | |
165 | cfContainer->SetVarTitle(idim, histoTitle.Data()); | |
166 | cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]); | |
9728bcfd | 167 | } |
a07db2fb | 168 | |
169 | TString stepTitle[kNsteps] = {"reconstructed", "generated"}; | |
589e3f71 | 170 | |
a07db2fb | 171 | TAxis* currAxis = 0x0; |
589e3f71 | 172 | for (Int_t istep=0; istep<kNsteps; istep++){ |
a07db2fb | 173 | cfContainer->SetStepTitle(istep, stepTitle[istep].Data()); |
174 | AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep); | |
175 | ||
176 | currAxis = gridSparse->GetAxis(kHvarMotherType); | |
177 | for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) { | |
178 | currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName()); | |
9728bcfd | 179 | } |
180 | } | |
4ab8d5a6 | 181 | |
a07db2fb | 182 | AddObjectToCollection(cfContainer, kTrackContainer); |
183 | ||
c6e0141f | 184 | fMuonTrackCuts->Print("mask"); |
aad6618e | 185 | } |
186 | ||
187 | //________________________________________________________________________ | |
a07db2fb | 188 | void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality) |
aad6618e | 189 | { |
190 | // | |
a07db2fb | 191 | /// Fill output objects |
aad6618e | 192 | // |
193 | ||
12e33589 | 194 | if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return; |
b201705a | 195 | |
12e33589 | 196 | Double_t ipVz = GetVertexSPD()->GetZ(); |
a07db2fb | 197 | Double_t ipVzMC = 0; |
198 | if ( IsMC() ) { | |
199 | if ( fMCEvent ) ipVzMC = fMCEvent->GetPrimaryVertex()->GetZ(); | |
7e7619d1 | 200 | else if ( fAODEvent ) { |
a07db2fb | 201 | AliAODMCHeader* aodMCHeader = (AliAODMCHeader *)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName()); |
202 | if ( aodMCHeader ) ipVzMC = aodMCHeader->GetVtxZ(); | |
55065f3f | 203 | } |
589e3f71 | 204 | } |
a07db2fb | 205 | |
206 | for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) { | |
207 | TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString(); | |
208 | ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz); | |
9728bcfd | 209 | } |
aad6618e | 210 | |
a07db2fb | 211 | 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 |
212 | if ( isPileupFromSPD ) return; | |
213 | ||
214 | Double_t containerInput[kNvars]; | |
215 | AliVParticle* track = 0x0; | |
93d6def1 | 216 | |
a07db2fb | 217 | for ( Int_t istep = 0; istep<2; ++istep ) { |
218 | Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks(); | |
219 | for (Int_t itrack = 0; itrack < nTracks; itrack++) { | |
220 | track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack); | |
221 | ||
c6e0141f | 222 | Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 ); |
a07db2fb | 223 | if ( ! isSelected ) continue; |
224 | ||
c6e0141f | 225 | // In W simulations with Pythia, sometimes muon is stored twice. |
226 | // Remove muon in case it has another muon as daugther | |
227 | if ( istep == kStepGeneratedMC ) { | |
228 | Int_t firstDaughter = GetDaughterIndex(track, 0); | |
229 | if ( firstDaughter >= 0 ) { | |
230 | Bool_t hasMuonDaughter = kFALSE; | |
231 | Int_t lastDaughter = GetDaughterIndex(track, 1); | |
232 | for ( Int_t idaugh=firstDaughter; idaugh<=lastDaughter; idaugh++ ) { | |
233 | AliVParticle* currTrack = GetMCTrack(idaugh); | |
234 | if ( currTrack->PdgCode() == track->PdgCode() ) { | |
235 | hasMuonDaughter = kTRUE; | |
236 | break; | |
237 | } | |
238 | } | |
239 | if ( hasMuonDaughter ) { | |
240 | AliDebug(1, Form("Current muon (%i) has muon daughter: rejecting it", itrack)); | |
241 | continue; | |
242 | } | |
243 | } | |
244 | } | |
245 | ||
a07db2fb | 246 | Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track); |
247 | ||
248 | Double_t thetaAbsEndDeg = 0; | |
249 | if ( istep == kStepReconstructed ) { | |
250 | Double_t rAbsEnd = ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd(); | |
251 | thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg(); | |
93d6def1 | 252 | } |
a07db2fb | 253 | else { |
254 | thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg(); | |
93d6def1 | 255 | } |
a07db2fb | 256 | Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310; |
257 | ||
258 | containerInput[kHvarPt] = track->Pt(); | |
259 | containerInput[kHvarEta] = track->Eta(); | |
260 | containerInput[kHvarPhi] = track->Phi(); | |
261 | containerInput[kHvarVz] = ( istep == kStepReconstructed ) ? ipVz : ipVzMC; | |
262 | containerInput[kHvarCharge] = track->Charge()/3.; | |
263 | containerInput[kHvarThetaAbs] = (Double_t)thetaAbsBin; | |
264 | containerInput[kHvarMotherType] = (Double_t)trackSrc; | |
265 | ||
266 | for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) { | |
267 | TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString(); | |
268 | if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue; | |
269 | ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep); | |
270 | } // loop on selected trigger classes | |
271 | } // loop on tracks | |
272 | } // loop on container steps | |
589e3f71 | 273 | } |
662e37fe | 274 | |
662e37fe | 275 | |
9728bcfd | 276 | //________________________________________________________________________ |
a07db2fb | 277 | void AliAnalysisTaskSingleMu::Terminate(Option_t *) { |
589e3f71 | 278 | // |
a07db2fb | 279 | /// Draw some histograms at the end. |
589e3f71 | 280 | // |
589e3f71 | 281 | |
a07db2fb | 282 | AliVAnalysisMuon::Terminate(""); |
b201705a | 283 | |
a07db2fb | 284 | if ( ! fMergeableCollection ) return; |
285 | ||
286 | TString physSel = fTerminateOptions->At(0)->GetName(); | |
287 | TString trigClassName = fTerminateOptions->At(1)->GetName(); | |
288 | TString centralityRange = fTerminateOptions->At(2)->GetName(); | |
289 | TString furtherOpt = fTerminateOptions->At(3)->GetName(); | |
12e33589 | 290 | |
291 | TString minBiasTrig = ""; | |
292 | TObjArray* optArr = furtherOpt.Tokenize(" "); | |
293 | TString currName = ""; | |
294 | for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) { | |
295 | currName = optArr->At(iopt)->GetName(); | |
df1e3f7a | 296 | if ( currName.Contains("-B-") || currName.Contains("ANY") ) minBiasTrig = currName; |
12e33589 | 297 | } |
298 | delete optArr; | |
299 | ||
a07db2fb | 300 | furtherOpt.ToUpper(); |
df1e3f7a | 301 | Bool_t plotChargeAsymmetry = furtherOpt.Contains("ASYM"); |
a07db2fb | 302 | |
303 | AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") ); | |
304 | if ( ! cfContainer ) return; | |
c6e0141f | 305 | |
306 | AliCFEffGrid* effSparse = new AliCFEffGrid(Form("eff%s", cfContainer->GetName()),Form("Efficiency %s", cfContainer->GetTitle()),*cfContainer); | |
307 | effSparse->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC); | |
308 | ||
309 | AliCFGridSparse* gridSparseArray[3] = {effSparse->GetNum(), effSparse->GetDen(), effSparse}; | |
310 | TString gridSparseName[3] = {cfContainer->GetStepTitle(kStepReconstructed), cfContainer->GetStepTitle(kStepGeneratedMC), "Efficiency"}; | |
a07db2fb | 311 | |
c6e0141f | 312 | Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange}; |
a07db2fb | 313 | // TString allSrcNames = ""; |
314 | // for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) { | |
315 | // if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" "); | |
316 | // allSrcNames += fSrcKeys->At(isrc)->GetName(); | |
317 | // } | |
318 | ||
319 | TCanvas* can = 0x0; | |
320 | Int_t xshift = 100; | |
321 | Int_t yshift = 100; | |
322 | Int_t igroup1 = -1; | |
323 | Int_t igroup2 = 0; | |
324 | ||
325 | Bool_t isMC = furtherOpt.Contains("MC"); | |
326 | Int_t firstSrc = ( isMC ) ? 0 : kUnidentified; | |
327 | Int_t lastSrc = ( isMC ) ? kNtrackSources - 1 : kUnidentified; | |
328 | if ( ! isMC ) srcColors[kUnidentified] = 1; | |
329 | ||
12e33589 | 330 | TString histoName = "", histoPattern = "", drawOpt = ""; |
a07db2fb | 331 | //////////////// |
332 | // Kinematics // | |
333 | //////////////// | |
c6e0141f | 334 | TCanvas* canKine[3] = {0x0, 0x0, 0x0}; |
335 | TLegend* legKine[3] = {0x0, 0x0, 0x0}; | |
336 | for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) { | |
df1e3f7a | 337 | for ( Int_t icharge=0; icharge<2; ++icharge ) { |
c6e0141f | 338 | for ( Int_t igrid=0; igrid<3; ++igrid ) { |
339 | if ( gridSparseArray[igrid]->GetEntries() == 0. ) break; | |
340 | if ( gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) { | |
341 | SetSparseRange(gridSparseArray[igrid], kHvarEta, "", -3.999, -2.501); | |
342 | SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc+1, isrc+1, "USEBIN"); | |
343 | SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge+1, icharge+1, "USEBIN"); | |
344 | } | |
345 | if ( ! canKine[igrid] ) { | |
346 | igroup1++; | |
347 | igroup2 = 0; | |
348 | currName = Form("%s_proj_%s", GetName(), gridSparseName[igrid].Data()); | |
349 | canKine[igrid] = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); | |
350 | canKine[igrid]->Divide(2,2); | |
351 | legKine[igrid] = new TLegend(0.6, 0.6, 0.8, 0.8); | |
352 | igroup2++; | |
353 | } | |
354 | for ( Int_t iproj=0; iproj<4; ++iproj ) { | |
355 | canKine[igrid]->cd(iproj+1); | |
356 | if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy(); | |
357 | TH1* projHisto = gridSparseArray[igrid]->Project(iproj); | |
df1e3f7a | 358 | projHisto->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), fSrcKeys->At(isrc)->GetName(), fChargeKeys->At(icharge)->GetName())); |
a07db2fb | 359 | if ( projHisto->GetEntries() == 0 ) continue; |
360 | Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 ); | |
c6e0141f | 361 | drawOpt = isFirst ? "e" : "esames"; |
a07db2fb | 362 | //if ( isrc == kUnidentified && ! drawOpt.Contains("same") ) isMC = kFALSE; |
363 | //if ( ! isMC ) srcColors[kUnidentified] = 1; | |
364 | projHisto->SetLineColor(srcColors[isrc]); | |
365 | projHisto->SetMarkerColor(srcColors[isrc]); | |
366 | projHisto->SetMarkerStyle(20+4*icharge); | |
367 | projHisto->Draw(drawOpt.Data()); | |
c6e0141f | 368 | gPad->Update(); |
369 | TPaveStats* paveStats = (TPaveStats*)projHisto->FindObject("stats"); | |
370 | if ( paveStats ) paveStats->SetTextColor(srcColors[isrc]); | |
a07db2fb | 371 | if ( iproj == 0 ) { |
372 | TString legEntry = fChargeKeys->At(icharge)->GetName(); | |
373 | if ( isMC ) legEntry += Form(" %s", fSrcKeys->At(isrc)->GetName()); | |
c6e0141f | 374 | legKine[igrid]->AddEntry(projHisto,legEntry.Data(), "lp"); |
a07db2fb | 375 | } |
c6e0141f | 376 | } // loop on grid sparse |
377 | } // loop on projections | |
378 | } // loop on mu charge | |
379 | } // loop on track sources | |
380 | ||
381 | ||
382 | for ( Int_t igrid=0; igrid<3; igrid++ ) { | |
383 | if ( ! canKine[igrid] ) continue; | |
384 | canKine[igrid]->cd(1); | |
385 | legKine[igrid]->Draw("same"); | |
386 | if ( gridSparseArray[igrid]->IsA() == AliCFEffGrid::Class() ) continue; | |
387 | SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", 1, gridSparseArray[igrid]->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range | |
a07db2fb | 388 | } // loop on container steps |
df1e3f7a | 389 | |
390 | // Plot charge asymmetry or mu+/mu- | |
391 | TString basePlotName = plotChargeAsymmetry ? "ChargeAsymmetry" : "ChargeRatio"; | |
392 | for ( Int_t igrid=0; igrid<2; igrid++ ) { | |
393 | if ( ! canKine[igrid] ) continue; | |
394 | TList* padList = canKine[igrid]->GetListOfPrimitives(); | |
395 | currName = canKine[igrid]->GetName(); | |
396 | currName.Append(Form("_%s", basePlotName.Data())); | |
397 | can = new TCanvas(currName.Data(),currName.Data(),canKine[igrid]->GetWindowTopX(),igroup2*yshift,600,600); | |
398 | can->Divide(2,2); | |
399 | for ( Int_t ipad=0; ipad<padList->GetEntries(); ipad++ ) { | |
400 | TPad* pad = dynamic_cast<TPad*> (padList->At(ipad)); | |
401 | if ( ! pad ) continue; | |
402 | TList* histoList = pad->GetListOfPrimitives(); | |
403 | can->cd(ipad+1); | |
404 | for ( Int_t iobj=0; iobj<histoList->GetEntries(); iobj++ ) { | |
405 | currName = histoList->At(iobj)->GetName(); | |
406 | if ( ! histoList->At(iobj)->InheritsFrom(TH1::Class()) || ! currName.Contains(fChargeKeys->At(1)->GetName()) ) continue; | |
407 | histoName = currName; | |
408 | histoName.ReplaceAll(fChargeKeys->At(1)->GetName(),""); | |
409 | histoName.Append(Form("_%s", basePlotName.Data())); | |
410 | currName.ReplaceAll(fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName()); | |
411 | TH1* auxHisto = dynamic_cast<TH1*> (histoList->FindObject(currName.Data())); | |
412 | if ( ! auxHisto ) continue; | |
413 | TH1* histo = static_cast<TH1*> (histoList->At(iobj)->Clone(histoName.Data())); | |
414 | if ( plotChargeAsymmetry ) { | |
415 | histo->Add(auxHisto, -1.); | |
416 | // h2 + h1 = 2xh2 + (h1-h2) | |
417 | auxHisto->Add(auxHisto, histo, 2.); | |
418 | } | |
419 | histo->Divide(auxHisto); | |
420 | histo->SetMarkerStyle(20); | |
421 | 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()); | |
422 | axisTitle.ReplaceAll("MuPlus","#mu^{+}"); | |
423 | axisTitle.ReplaceAll("MuMinus","#mu^{-}"); | |
424 | histo->GetYaxis()->SetTitle(axisTitle.Data()); | |
425 | histo->SetStats(kFALSE); | |
426 | drawOpt = ( gPad->GetListOfPrimitives()->GetEntries() == 0 ) ? "e" : "esames"; | |
427 | histo->Draw(drawOpt.Data()); | |
428 | } // loop on histos | |
429 | gPad->Update(); | |
430 | } // loop on pads | |
431 | } // loop on container steps | |
a07db2fb | 432 | |
433 | ||
434 | ////////////////////// | |
435 | // Event statistics // | |
436 | ////////////////////// | |
437 | printf("\nTotal analyzed events:\n"); | |
438 | TString evtSel = Form("trigger:%s", trigClassName.Data()); | |
439 | fEventCounters->PrintSum(evtSel.Data()); | |
440 | printf("Physics selected analyzed events:\n"); | |
441 | evtSel = Form("trigger:%s/selected:yes", trigClassName.Data()); | |
442 | fEventCounters->PrintSum(evtSel.Data()); | |
443 | ||
444 | TString countPhysSel = "any"; | |
445 | if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes"; | |
446 | else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no"; | |
447 | countPhysSel.Prepend("selected:"); | |
448 | printf("Analyzed events vs. centrality:\n"); | |
449 | evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data()); | |
450 | fEventCounters->Print("centrality",evtSel.Data(),kTRUE); | |
451 | ||
452 | ||
453 | /////////////////// | |
454 | // Vertex method // | |
455 | /////////////////// | |
456 | if ( ! furtherOpt.Contains("VERTEX") ) return; | |
df1e3f7a | 457 | Int_t firstMother = ( isMC ) ? 0 : kUnidentified; |
458 | Int_t lastMother = ( isMC ) ? kNtrackSources - 1 : kUnidentified; | |
a07db2fb | 459 | igroup1++; |
12e33589 | 460 | TH1* eventVertex = (TH1*)GetSum(physSel, minBiasTrig, centralityRange, "hIpVtx"); |
a07db2fb | 461 | if ( ! eventVertex ) return; |
462 | Double_t minZ = -9.99, maxZ = 9.99; | |
463 | Double_t meanZ = 0., sigmaZ = 4.; | |
464 | Double_t nSigma = 2.; | |
12e33589 | 465 | TString fitOpt = "R0S"; |
a07db2fb | 466 | Bool_t fixFitRange = kFALSE; |
467 | TString fitFormula = Form("[0]+[1]*(x+[2])"); | |
468 | ||
469 | // Get vertex shape | |
470 | if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2(); | |
471 | Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow | |
472 | printf("Event vertex integral %.0f\n\n", eventVtxIntegral); | |
473 | if ( eventVtxIntegral <= 0. ) return; | |
474 | eventVertex->Scale(1./eventVtxIntegral); | |
475 | printf("\nFit MB vertex\n"); | |
476 | eventVertex->Fit("gaus",fitOpt.Data(),"",minZ,maxZ); | |
477 | TF1* vtxFit = (TF1*)eventVertex->GetListOfFunctions()->FindObject("gaus"); | |
478 | currName = "vtxIntegrated"; | |
479 | can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); | |
480 | can->SetLogy(); | |
481 | eventVertex->Draw(); | |
482 | vtxFit->Draw("same"); | |
b201705a | 483 | |
a07db2fb | 484 | |
485 | enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos}; | |
486 | TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"}; | |
487 | TArrayI sumMothers[kNrecoHistos]; | |
488 | sumMothers[kRecoHF].Set(0); | |
489 | sumMothers[kRecoBkg].Set(0); | |
490 | sumMothers[kInputHF].Set(3); | |
491 | sumMothers[kInputHF][0] = kCharmMu; | |
492 | sumMothers[kInputHF][1] = kBeautyMu; | |
493 | sumMothers[kInputHF][2] = kQuarkoniumMu; | |
494 | sumMothers[kInputDecay].Set(1); | |
495 | sumMothers[kInputDecay][0] = kDecayMu; | |
496 | sumMothers[kRecoAll].Set(kNtrackSources); | |
497 | for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) { | |
498 | sumMothers[kRecoAll][isrc] = isrc; | |
b201705a | 499 | } |
a07db2fb | 500 | |
501 | meanZ = vtxFit->GetParameter(1); | |
502 | sigmaZ = vtxFit->GetParameter(2); | |
503 | ||
504 | Double_t minZfit = ( fixFitRange ) ? minZ : meanZ - nSigma*sigmaZ; | |
505 | Double_t maxZfit = ( fixFitRange ) ? maxZ : meanZ + nSigma*sigmaZ; | |
506 | ||
507 | TF1* fitFunc = new TF1("fitFunc", fitFormula.Data(), minZ, maxZ); | |
508 | fitFunc->SetLineColor(2); | |
509 | fitFunc->SetParNames("Line norm", "Line slope", "Free path"); | |
510 | const Double_t kFreePath = 153.; // 150.; // 130.; // cm | |
511 | //fitFunc->SetParameters(0.,1.); | |
512 | fitFunc->FixParameter(2, kFreePath); | |
513 | ||
514 | AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed); | |
515 | TAxis* ptAxis = gridSparse->GetAxis(kHvarPt); | |
516 | ||
517 | Double_t slope = 0.; | |
518 | Double_t limitNorm = 0., limitSlope = 0.; | |
12e33589 | 519 | Int_t firstPtBin = 0, lastPtBin = 0; |
520 | ||
521 | gStyle->SetOptFit(1111); | |
a07db2fb | 522 | |
523 | for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) { | |
524 | igroup2++; | |
525 | SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN"); | |
526 | SetSparseRange(gridSparse, kHvarPt, "", 1, ptAxis->GetNbins(), "USEBIN"); | |
527 | TH1* recoHisto[kNrecoHistos]; | |
528 | for ( Int_t ireco=0; ireco<kNrecoHistos; ++ireco ) { | |
529 | recoHisto[ireco] = gridSparse->Project(kHvarPt); | |
530 | histoName = Form("%sMuon_%s", baseRecoName[ireco].Data(), fThetaAbsKeys->At(itheta)->GetName()); | |
531 | recoHisto[ireco]->SetName(histoName.Data()); | |
532 | recoHisto[ireco]->SetTitle(histoName.Data()); | |
533 | recoHisto[ireco]->Reset(); | |
534 | recoHisto[ireco]->Sumw2(); | |
535 | for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) { | |
536 | SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc]+1, sumMothers[ireco][isrc]+1, "USEBIN"); | |
537 | TH1* auxHisto = gridSparse->Project(kHvarPt); | |
538 | recoHisto[ireco]->Add(auxHisto); | |
539 | delete auxHisto; | |
540 | } | |
b201705a | 541 | } |
a07db2fb | 542 | SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN"); |
543 | Int_t currDraw = 0; | |
544 | ||
545 | for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) { | |
546 | firstPtBin = ibinpt; | |
547 | lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt; | |
548 | SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN"); | |
549 | TH1* histo = gridSparse->Project(kHvarVz); | |
550 | histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt)); | |
551 | if ( histo->GetEntries() < 100. ) break; | |
552 | printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries()); | |
553 | histo->Divide(eventVertex); | |
554 | Double_t norm = histo->GetBinContent(histo->FindBin(0.)); | |
555 | histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)"); | |
556 | slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) - | |
557 | histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ ); | |
558 | ||
559 | if ( slope < 0. ) slope = norm/kFreePath; | |
560 | ||
561 | // Try to fit twice: it fit fails the first time | |
562 | // set some limits on parameters | |
563 | for ( Int_t itry=0; itry<2; itry++ ) { | |
564 | fitFunc->SetParameter(0, norm); | |
565 | fitFunc->SetParameter(1, slope); | |
566 | if ( itry > 0 ) { | |
567 | limitNorm = 2.*histo->Integral(); | |
568 | limitSlope = 2.*histo->Integral()/kFreePath; | |
569 | //fitFunc->SetParLimits(0, 0., limitNorm); // REMEMBER TO CHECK | |
570 | fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK | |
571 | printf("Norm 0. < %f < %f slope 0. < %f < %f\n", norm, limitNorm, slope, limitSlope); | |
572 | } | |
12e33589 | 573 | TFitResultPtr fitRes = histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit); |
a07db2fb | 574 | |
12e33589 | 575 | // if ( gMinuit->fCstatu.Contains("CONVERGED") && |
576 | if ( ((Int_t)fitRes) == 0 && | |
a07db2fb | 577 | fitFunc->GetParameter(0) > 0. && |
578 | fitFunc->GetParameter(1) > 0. ) | |
579 | break; | |
7e7619d1 | 580 | else if ( furtherOpt.Contains("REFIT") ) printf("Re-fit with limits\n"); |
a07db2fb | 581 | else { |
7e7619d1 | 582 | printf("Warning: fit problems !!!\n"); |
583 | break; | |
a07db2fb | 584 | } |
585 | } | |
586 | ||
587 | Double_t p0 = fitFunc->GetParameter(0); | |
588 | Double_t p0err = fitFunc->GetParError(0); | |
589 | Double_t p1 = fitFunc->GetParameter(1); | |
590 | Double_t p1err = fitFunc->GetParError(1); | |
591 | ||
592 | Double_t nullVz = ( p1 != 0. ) ? -p0/p1 : 0.; | |
593 | Double_t nullVzErr = ( p0 != 0. && p1 != 0. ) ? TMath::Abs(nullVz) * TMath::Sqrt(p0err*p0err/(p0*p0) + p1err*p1err/(p1*p1) ) : 0.; | |
594 | ||
595 | printf("Null value at %f +- %f\n", nullVz - kFreePath, nullVzErr); | |
596 | ||
597 | recoHisto[kRecoHF]->SetBinContent(ibinpt, p0); | |
598 | recoHisto[kRecoHF]->SetBinError(ibinpt, p0err); | |
599 | recoHisto[kRecoBkg]->SetBinContent(ibinpt, ( kFreePath + meanZ ) * p1); | |
600 | recoHisto[kRecoBkg]->SetBinError(ibinpt, ( kFreePath + meanZ ) * p1err); | |
601 | if ( currDraw%4 == 0 ){ | |
602 | currName = Form("vtx_%s_PtBin%i",fThetaAbsKeys->At(itheta)->GetName(), ibinpt); | |
603 | can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); | |
604 | can->Divide(2,2); | |
605 | } | |
606 | can->cd( currDraw%4 + 1 ); | |
607 | can->SetLogy(); | |
608 | histo->Draw(); | |
609 | fitFunc->DrawCopy("same"); | |
610 | currDraw++; | |
611 | } // loop on pt bins | |
612 | SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN"); | |
a07db2fb | 613 | currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName()); |
614 | can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600); | |
615 | TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8); | |
616 | drawOpt = "e"; | |
617 | for ( Int_t ireco=0; ireco<kNrecoHistos-1; ++ireco ) { | |
618 | if ( recoHisto[ireco]->GetEntries() == 0. ) continue; | |
619 | TH1* ratio = (TH1*)recoHisto[ireco]->Clone(Form("%s_ratio", recoHisto[ireco]->GetName())); | |
620 | ratio->Divide(recoHisto[kRecoAll]); | |
621 | ratio->SetLineColor(srcColors[ireco]); | |
622 | ratio->SetMarkerColor(srcColors[ireco]); | |
623 | ratio->SetMarkerStyle(20+ireco); | |
624 | ratio->GetYaxis()->SetTitle("fraction of total"); | |
625 | ratio->Draw(drawOpt.Data()); | |
626 | leg->AddEntry(ratio,baseRecoName[ireco].Data(), "lp"); | |
627 | drawOpt = "esame"; | |
b201705a | 628 | } |
a07db2fb | 629 | leg->Draw("same"); |
630 | } // loop on theta abs | |
589e3f71 | 631 | } |