Fix problem in terminate. Add possibility to select the minimum number of vertex...
[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisTaskSingleMu.cxx
CommitLineData
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
77ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
78/// \endcond
aad6618e 79
9728bcfd 80
aad6618e 81//________________________________________________________________________
a07db2fb 82AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() :
83 AliVAnalysisMuon(),
84 fThetaAbsKeys(0x0)
aad6618e 85{
a07db2fb 86 /// Default ctor.
aad6618e 87}
88
9728bcfd 89//________________________________________________________________________
a07db2fb 90AliAnalysisTaskSingleMu::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//________________________________________________________________________
103AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
589e3f71 104{
105 //
a07db2fb 106 /// Destructor
589e3f71 107 //
108
a07db2fb 109 delete fThetaAbsKeys;
aad6618e 110}
111
112//___________________________________________________________________________
a07db2fb 113void 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 188void 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 277void 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}