]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/muon/AliAnalysisTaskSingleMu.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates (F. Bellini)
[u/mrichter/AliRoot.git] / PWG3 / 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
662e37fe 16//-----------------------------------------------------------------------------
17/// \class AliAnalysisTaskSingleMu
18/// Analysis task for single muons in the spectrometer.
9728bcfd 19/// The output is a list of histograms.
b201705a 20/// The macro class can run on AODs or ESDs.
21/// In the latter case a flag can be activated to produce a tree as output.
9728bcfd 22/// If Monte Carlo information is present, some basics checks are performed.
662e37fe 23///
24/// \author Diego Stocco
25//-----------------------------------------------------------------------------
26
aad6618e 27//----------------------------------------------------------------------------
28// Implementation of the single muon analysis class
662e37fe 29// An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
aad6618e 30//----------------------------------------------------------------------------
31
32#define AliAnalysisTaskSingleMu_cxx
33
34// ROOT includes
aad6618e 35#include "TROOT.h"
9728bcfd 36#include "TH1.h"
37#include "TH2.h"
38#include "TH3.h"
39#include "TAxis.h"
40#include "TList.h"
662e37fe 41#include "TCanvas.h"
9728bcfd 42#include "TMath.h"
b201705a 43#include "TTree.h"
44#include "TTimeStamp.h"
589e3f71 45#include "TMap.h"
46#include "TObjString.h"
55065f3f 47#include "TObjArray.h"
589e3f71 48#include "TIterator.h"
49#include "TParameter.h"
50#include "TMCProcess.h"
aad6618e 51
52// STEER includes
53#include "AliLog.h"
54
589e3f71 55#include "AliAODInputHandler.h"
aad6618e 56#include "AliAODEvent.h"
57#include "AliAODTrack.h"
58#include "AliAODVertex.h"
aad6618e 59
9728bcfd 60#include "AliMCParticle.h"
61#include "AliStack.h"
62
b201705a 63#include "AliESDInputHandler.h"
64#include "AliESDEvent.h"
65#include "AliESDMuonTrack.h"
9728bcfd 66
589e3f71 67#include "AliMultiplicity.h"
68
69// ANALYSIS includes
70#include "AliAnalysisManager.h"
9728bcfd 71#include "AliAnalysisTaskSE.h"
589e3f71 72
73// CORRFW includes
74#include "AliCFManager.h"
75#include "AliCFContainer.h"
76#include "AliCFGridSparse.h"
77#include "AliCFTrackKineCuts.h"
78#include "AliCFParticleGenCuts.h"
79#include "AliCFEventRecCuts.h"
80
aad6618e 81#include "AliAnalysisTaskSingleMu.h"
82
662e37fe 83/// \cond CLASSIMP
84ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
85/// \endcond
aad6618e 86
9728bcfd 87
aad6618e 88//________________________________________________________________________
589e3f71 89AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, Int_t fillTreeScaleDown, Bool_t keepAll) :
9728bcfd 90 AliAnalysisTaskSE(name),
589e3f71 91 fFillTreeScaleDown(fillTreeScaleDown),
b201705a 92 fKeepAll(keepAll),
589e3f71 93 fkNvtxContribCut(1),
94 fCFManager(0),
9728bcfd 95 fHistoList(0),
b201705a 96 fHistoListMC(0),
589e3f71 97 fHistoListQA(0),
b201705a 98 fTreeSingleMu(0),
99 fTreeSingleMuMC(0),
100 fVarFloat(0),
101 fVarInt(0),
102 fVarChar(0),
103 fVarUInt(0),
104 fVarFloatMC(0),
589e3f71 105 fVarIntMC(0),
55065f3f 106 fAuxObjects(new TMap()),
107 fTriggerClasses(0x0)
aad6618e 108{
109 //
110 /// Constructor.
111 //
589e3f71 112 if ( fFillTreeScaleDown <= 0 )
b201705a 113 fKeepAll = kFALSE;
114
589e3f71 115 DefineOutput(1, AliCFContainer::Class());
9728bcfd 116 DefineOutput(2, TList::Class());
589e3f71 117 DefineOutput(3, TList::Class());
118 DefineOutput(4, TList::Class());
119
120 if ( fFillTreeScaleDown > 0 )
121 DefineOutput(5, TTree::Class());
122
55065f3f 123 fAuxObjects->SetOwner();
b201705a 124
55065f3f 125 SetTriggerClasses();
aad6618e 126}
127
9728bcfd 128
129//________________________________________________________________________
130AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
aad6618e 131{
132 //
9728bcfd 133 /// Destructor
aad6618e 134 //
589e3f71 135
136 delete fCFManager->GetParticleContainer(); // The container is not deleted by framework
137 delete fCFManager;
b201705a 138 delete fHistoList;
139 delete fHistoListMC;
589e3f71 140 delete fHistoListQA;
b201705a 141 delete fTreeSingleMu;
142 delete fTreeSingleMuMC;
143 delete fVarFloat;
144 delete fVarInt;
145 delete [] fVarChar;
146 delete fVarUInt;
147 delete fVarFloatMC;
148 delete fVarIntMC;
55065f3f 149 delete fAuxObjects;
150 delete fTriggerClasses;
9728bcfd 151}
aad6618e 152
9728bcfd 153
154//___________________________________________________________________________
b201705a 155void AliAnalysisTaskSingleMu::NotifyRun()
9728bcfd 156{
157 //
158 /// Use the event handler information to correctly fill the analysis flags:
159 /// - check if Monte Carlo information is present
160 //
161
589e3f71 162 if ( fMCEvent ) {
9728bcfd 163 AliInfo("Monte Carlo information is present");
164 }
165 else {
166 AliInfo("No Monte Carlo information in run");
aad6618e 167 }
589e3f71 168
169 Int_t histoIndex = -1;
170
171 TString histoName = Form("hVtx%u",InputEvent()->GetRunNumber());
55065f3f 172 if ( ! fAuxObjects->GetValue(histoName.Data()) ) {
589e3f71 173 histoIndex = GetHistoIndex(kHistoEventVz);
174 AliInfo(Form("Creating histogram %s", histoName.Data()));
175 TH1F* histo1D = (TH1F*)fHistoList->At(histoIndex)->Clone(histoName.Data());
176 histo1D->Reset();
177 histo1D->SetNormFactor(0.);
178 histo1D->Sumw2();
55065f3f 179 fAuxObjects->Add(new TObjString(histoName.Data()), histo1D);
589e3f71 180 }
181
182 /*
183 histoIndex = GetHistoIndex(kHistoNmuonsPerRun);
184 TH2F* histoMultPerRun = (TH2F*)fHistoList->At(histoIndex);
185 histoMultPerRun->GetXaxis()->SetBinLabel(1,Form("%u",InputEvent()->GetRunNumber()));
186 */
187}
188
189
190//___________________________________________________________________________
191void AliAnalysisTaskSingleMu::FinishTaskOutput()
192{
193 //
194 /// Perform histogram normalization after the last analyzed event
195 /// but before merging
196 //
197
198 // Normalize the vertex distribution histogram
199 // to the total number of analyzed muons
200 Int_t histoIndex = GetHistoIndex(kHistoEventVzNorm);
201 TH1F* histoVtxAll = (TH1F*)fHistoList->At(histoIndex);
55065f3f 202 TIter nextVal(fAuxObjects);
589e3f71 203 TObjString* objString = 0x0;
204 while( ( objString = (TObjString*)nextVal() ) ){
55065f3f 205 if ( ! objString->GetString().Contains("hVtx") ) continue;
206 TH1F* currHisto = (TH1F*)fAuxObjects->GetValue(objString->String().Data());
589e3f71 207 Double_t mbEventsPerRun = currHisto->Integral();
208 if ( mbEventsPerRun == 0. )
209 continue;
210 Double_t muEventsPerRun = currHisto->GetNormFactor();
211 if ( muEventsPerRun == 0. )
212 continue;
213 AliInfo(Form("Add vertex histogram %s. Events: MB %i Mu %i",currHisto->GetName(), (Int_t)mbEventsPerRun, (Int_t)muEventsPerRun));
214 histoVtxAll->Add(currHisto, muEventsPerRun/mbEventsPerRun);
215 }
216
217 // Set the correct run limits for the histograms
218 // vs run number to cope with the merging
219 Int_t indexPerRun[2] = {kHistoNeventsPerRun, kHistoNmuonsPerRun};
220 for ( Int_t ihisto=0; ihisto<2; ihisto++) {
221 histoIndex = GetHistoIndex(indexPerRun[ihisto]);
222 TH2F* histo2D = (TH2F*)fHistoList->At(histoIndex);
223 Double_t minX = 1e10, maxX = -1e10;
224 for (Int_t ibin=1; ibin<=histo2D->GetXaxis()->GetNbins(); ibin++){
225 TString runNum = histo2D->GetXaxis()->GetBinLabel(ibin);
226 minX = TMath::Min(runNum.Atof()-0.5, minX);
227 maxX = TMath::Max(runNum.Atof()+0.5, maxX);
228 }
229 histo2D->GetXaxis()->SetLimits(minX, maxX);
230 AliInfo(Form("Histogram %s run limits (%f, %f)",histo2D->GetName(), minX, maxX));
231 }
aad6618e 232}
233
9728bcfd 234
aad6618e 235//___________________________________________________________________________
9728bcfd 236void AliAnalysisTaskSingleMu::UserCreateOutputObjects()
aad6618e 237{
238 //
239 /// Create output histograms
240 //
9728bcfd 241 AliInfo(Form(" CreateOutputObjects of task %s\n", GetName()));
aad6618e 242
589e3f71 243 if ( fFillTreeScaleDown > 0 ) OpenFile(1);
b201705a 244
9728bcfd 245 // initialize histogram lists
589e3f71 246 if ( ! fHistoList ) fHistoList = new TList();
247 if ( ! fHistoListMC ) fHistoListMC = new TList();
248 if ( ! fHistoListQA ) fHistoListQA = new TList();
b201705a 249
250 // Init variables
251 fVarFloat = new Float_t [kNvarFloat];
252 fVarInt = new Int_t [kNvarInt];
253 fVarChar = new Char_t *[kNvarChar];
254 fVarUInt = new UInt_t [kNvarUInt];
255 fVarFloatMC = new Float_t [kNvarFloatMC];
256 fVarIntMC = new Int_t [kNvarIntMC];
589e3f71 257
258 const Int_t charWidth[kNvarChar] = {255};
259 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
260 fVarChar[ivar] = new Char_t [charWidth[ivar]];
261 }
262
9728bcfd 263 Int_t nPtBins = 60;
589e3f71 264 Double_t ptMin = 0., ptMax = 30.;
9728bcfd 265 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
589e3f71 266
55065f3f 267 Int_t nRapidityBins = 25;
268 Double_t rapidityMin = -4.5, rapidityMax = -2.;
269 TString rapidityName("Rapidity"), rapidityTitle("y"), rapidityUnits("a.u.");
589e3f71 270
271 Int_t nPhiBins = 36;
272 Double_t phiMin = 0.; Double_t phiMax = 2*TMath::Pi();
273 TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
9728bcfd 274
589e3f71 275 Int_t nDcaBins = 30;
276 Double_t dcaMin = 0., dcaMax = 300.;
9728bcfd 277 TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");
662e37fe 278
589e3f71 279 Int_t nVzBins = 40;
280 Double_t vzMin = -20., vzMax = 20.;
9728bcfd 281 TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
282
589e3f71 283 Int_t nThetaAbsEndBins = 4;
284 Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 3.5;
285 TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");
286
287 Int_t nChargeBins = 2;
288 Double_t chargeMin = -2, chargeMax = 2.;
289 TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
290
291 Int_t nMatchTrigBins = 4;
292 Double_t matchTrigMin = -0.5, matchTrigMax = 3.5;
293 TString matchTrigName("MatchTrig"), matchTrigTitle("Trigger match"), matchTrigUnits("");
9728bcfd 294
55065f3f 295 Int_t nTrigClassBins = fTriggerClasses->GetEntries();
296 Double_t trigClassMin = 0.5, trigClassMax = (Double_t)nTrigClassBins + 0.5;
589e3f71 297 TString tricClassName("TrigClass"), trigClassTitle("Fired trigger class"), trigClassUnits("");
298
299 Int_t nGoodVtxBins = 3;
300 Double_t goodVtxMin = -0.5, goodVtxMax = 2.5;
301 TString goodVtxName("GoodVtx"), goodVtxTitle("Vertex flags"), goodVtxUnits("");
302
303 Int_t nMotherTypeBins = kNtrackSources;
304 Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
305 TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
9728bcfd 306
307 TString trigName[kNtrigCuts];
308 trigName[kNoMatchTrig] = "NoMatch";
309 trigName[kAllPtTrig] = "AllPt";
310 trigName[kLowPtTrig] = "LowPt";
311 trigName[kHighPtTrig] = "HighPt";
312
313 TString srcName[kNtrackSources];
314 srcName[kCharmMu] = "Charm";
315 srcName[kBeautyMu] = "Beauty";
316 srcName[kPrimaryMu] = "Decay";
317 srcName[kSecondaryMu] = "Secondary";
318 srcName[kRecoHadron] = "Hadrons";
319 srcName[kUnknownPart] = "Unknown";
320
321 TH1F* histo1D = 0x0;
322 TH2F* histo2D = 0x0;
323 TString histoName, histoTitle;
324 Int_t histoIndex = 0;
325
589e3f71 326 // Multi-dimensional histo
55065f3f 327 Int_t nbins[kNvars] = {nPtBins, nRapidityBins, nPhiBins, nDcaBins, nVzBins, nThetaAbsEndBins, nChargeBins, nMatchTrigBins, nTrigClassBins, nGoodVtxBins, nMotherTypeBins};
328 Double_t xmin[kNvars] = {ptMin, rapidityMin, phiMin, dcaMin, vzMin, thetaAbsEndMin, chargeMin, matchTrigMin, trigClassMin, goodVtxMin, motherTypeMin};
329 Double_t xmax[kNvars] = {ptMax, rapidityMax, phiMax, dcaMax, vzMax, thetaAbsEndMax, chargeMax, matchTrigMax, trigClassMax, goodVtxMax, motherTypeMax};
330 TString axisTitle[kNvars] = {ptTitle, rapidityTitle, phiTitle, dcaTitle, vzTitle, thetaAbsEndTitle, chargeTitle, matchTrigTitle, trigClassTitle, goodVtxTitle, motherTypeTitle};
331 TString axisUnits[kNvars] = {ptUnits, rapidityUnits, phiUnits, dcaUnits, vzUnits, thetaAbsEndUnits, chargeUnits, matchTrigUnits, trigClassUnits, goodVtxUnits, motherTypeUnits};
662e37fe 332
55065f3f 333 TString stepTitle[kNsteps] = {"reconstructed", "in acceptance", "generated", "in acceptance (MC)"};
589e3f71 334
335 // Create CF container
336 AliCFContainer* container = new AliCFContainer("SingleMuonContainer","container for tracks",kNsteps,kNvars,nbins);
337
338 for ( Int_t idim = 0; idim<kNvars; idim++){
339 histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
340 histoTitle.ReplaceAll("()","");
341
342 container->SetVarTitle(idim, histoTitle.Data());
343 container->SetBinLimits(idim, xmin[idim], xmax[idim]);
9728bcfd 344 }
589e3f71 345
346 for (Int_t istep=0; istep<kNsteps; istep++){
347 container->SetStepTitle(istep, stepTitle[istep].Data());
348 AliCFGridSparse* gridSparse = container->GetGrid(istep);
349
350 SetAxisLabel(gridSparse->GetAxis(kHvarTrigClass));
351 TAxis* isGoodVtxAxis = gridSparse->GetAxis(kHvarIsGoodVtx);
352 isGoodVtxAxis->SetBinLabel(1,Form("No vertex (contrib<%i)",fkNvtxContribCut));
353 isGoodVtxAxis->SetBinLabel(2,"Good vertex");
354 isGoodVtxAxis->SetBinLabel(3,"Pileup");
355 TAxis* motherTypeAxis = gridSparse->GetAxis(kHvarMotherType);
356 for (Int_t ibin=0; ibin<kNtrackSources; ibin++){
357 motherTypeAxis->SetBinLabel(ibin+1,srcName[ibin]);
358 }
9728bcfd 359 }
589e3f71 360
361 // Create cuts
362
363 //Particle-Level cuts:
364 AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts();
365 mcGenCuts->SetNameTitle("mcGenCuts","MC particle generation cuts");
366 mcGenCuts->SetRequirePdgCode(13,kTRUE);
367 mcGenCuts->SetQAOn(fHistoListQA);
368
369 // MC kinematic cuts
370 AliCFTrackKineCuts *mcAccCuts = new AliCFTrackKineCuts();
371 mcAccCuts->SetNameTitle("mcAccCuts","MC-level acceptance cuts");
372 mcAccCuts->SetEtaRange(-4.,-2.5);
589e3f71 373 mcAccCuts->SetQAOn(fHistoListQA);
374
589e3f71 375 // Rec-Level kinematic cuts
376 AliCFTrackKineCuts *recAccCuts = new AliCFTrackKineCuts();
377 recAccCuts->SetNameTitle("recAccCuts","Reco-level acceptance cuts");
378 recAccCuts->SetEtaRange(-4.,-2.5);
379 recAccCuts->SetQAOn(fHistoListQA);
380
589e3f71 381 TObjArray* mcGenList = new TObjArray(0) ;
382 mcGenList->AddLast(mcGenCuts);
383
384 TObjArray* mcAccList = new TObjArray(0);
385 mcAccList->AddLast(mcAccCuts);
386
589e3f71 387 TObjArray* recAccList = new TObjArray(0);
388 recAccList->AddLast(recAccCuts);
389
589e3f71 390 // Create CF manager
391 fCFManager = new AliCFManager() ;
392 fCFManager->SetParticleContainer(container);
393
394 // Add cuts
395 // Dummy event container
396 Int_t dummyBins[1] = {1};
397 AliCFContainer* evtContainer = new AliCFContainer("dummyContainer","dummy contaier for events",1,1,dummyBins);
398 fCFManager->SetEventContainer(evtContainer);
399 fCFManager->SetEventCutsList(0,0x0);
400
401 // Init empty cuts (avoid warnings in framework)
402 for (Int_t istep=0; istep<kNsteps; istep++) {
403 fCFManager->SetParticleCutsList(istep,0x0);
662e37fe 404 }
589e3f71 405 fCFManager->SetParticleCutsList(kStepAcceptance,recAccList);
589e3f71 406 fCFManager->SetParticleCutsList(kStepGeneratedMC,mcGenList);
407 fCFManager->SetParticleCutsList(kStepAcceptanceMC,mcAccList);
662e37fe 408
b201705a 409 // Summary histos
589e3f71 410 histoName = "histoNeventsPerTrig";
411 histoTitle = "Number of events per trigger class";
412 histo2D = new TH2F(histoName.Data(), histoTitle.Data(), nTrigClassBins, trigClassMin, trigClassMax, 5, 0.5, 5.5);
413 histo2D->GetXaxis()->SetTitle("Fired class");
414 SetAxisLabel(histo2D->GetXaxis());
415 histo2D->GetYaxis()->SetBinLabel(1, "All events");
416 histo2D->GetYaxis()->SetBinLabel(2, "Good vtx events");
417 histo2D->GetYaxis()->SetBinLabel(3, "Pileup events");
418 histo2D->GetYaxis()->SetBinLabel(4, "All vertices");
419 histo2D->GetYaxis()->SetBinLabel(5, "Pileup vertices");
420 histoIndex = GetHistoIndex(kHistoNeventsPerTrig);
421 fHistoList->AddAt(histo2D, histoIndex);
422
b201705a 423 histoName = "histoMuonMultiplicity";
424 histoTitle = "Muon track multiplicity";
589e3f71 425 histo2D = new TH2F(histoName.Data(), histoTitle.Data(), 15, -0.5, 15-0.5,
426 nTrigClassBins, trigClassMin, trigClassMax);
427 histo2D->GetXaxis()->SetTitle("# of muons");
428 SetAxisLabel(histo2D->GetYaxis());
429 histoIndex = GetHistoIndex(kHistoMuonMultiplicity);
430 fHistoList->AddAt(histo2D, histoIndex);
431
432 histoName = "histoEventVz";
433 histoTitle = "All events IP Vz distribution";
434 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nVzBins, vzMin, vzMax);
435 histoTitle = Form("%s (%s)", vzTitle.Data(), vzUnits.Data());
436 histo1D->GetXaxis()->SetTitle(histoTitle.Data());
437 histoIndex = GetHistoIndex(kHistoEventVz);
438 fHistoList->AddAt(histo1D, histoIndex);
439
440 histoName = "histoEventVzNorm";
441 histoTitle = "All events IP Vz distribution normalized to CMUS1B";
442 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nVzBins, vzMin, vzMax);
443 histoTitle = Form("%s (%s)", vzTitle.Data(), vzUnits.Data());
444 histo1D->GetXaxis()->SetTitle(histoTitle.Data());
445 histoIndex = GetHistoIndex(kHistoEventVzNorm);
b201705a 446 fHistoList->AddAt(histo1D, histoIndex);
447
589e3f71 448 Int_t hRunIndex[2] = {kHistoNeventsPerRun, kHistoNmuonsPerRun};
449 TString hRunName[2] = {"histoNeventsPerRun", "histoNmuonsPerRun"};
450 TString hRunTitle[2] = {"Number of events per run", "Number of muons per run"};
451 for ( Int_t ihisto=0; ihisto<2; ihisto++){
452 histoName = hRunName[ihisto];
453 histoTitle = hRunTitle[ihisto];
454 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
455 1, 0., 0.,
456 nTrigClassBins, trigClassMin, trigClassMax);
457 histo2D->GetXaxis()->SetTitle("Run number");
458 SetAxisLabel(histo2D->GetYaxis());
459 histo2D->Sumw2();
460 histoIndex = hRunIndex[ihisto];
461 fHistoList->AddAt(histo2D, histoIndex);
462 }
463
464 // MC histos summary
465 histoName = "histoCheckVz";
466 histoTitle = "Check IP Vz distribution";
467 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
468 nVzBins, vzMin, vzMax,
469 nVzBins, vzMin, vzMax);
470 histoTitle = Form("%s (%s)", vzTitle.Data(), vzUnits.Data());
471 histo2D->GetXaxis()->SetTitle(histoTitle.Data());
472 histoTitle = Form("%s MC (%s)", vzTitle.Data(), vzUnits.Data());
473 histo2D->GetYaxis()->SetTitle(histoTitle.Data());
474 histoIndex = GetHistoIndex(kHistoCheckVzMC);
475 fHistoListMC->AddAt(histo2D, histoIndex);
476
9728bcfd 477 // MC histos
478 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
479 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
480 histoName = Form("%sResolution%sTrig%s", ptName.Data(), trigName[itrig].Data(), srcName[isrc].Data());
481 histoTitle = Form("%s resolution. Trigger: %s (%s)", ptTitle.Data(), trigName[itrig].Data(), srcName[isrc].Data());
482 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 100, -5., 5.);
483 histo1D->GetXaxis()->SetTitle(Form("%s^{reco} - %s^{MC} (%s)", ptTitle.Data(), ptTitle.Data(), ptUnits.Data()));
b201705a 484 histoIndex = GetHistoIndex(kHistoPtResolutionMC, itrig, isrc);
9728bcfd 485 fHistoListMC->AddAt(histo1D, histoIndex);
486 }
487 }
488
b201705a 489 // Trees
589e3f71 490 if ( fFillTreeScaleDown > 0 ){
b201705a 491 TString leavesFloat[kNvarFloat] = {"Px", "Py", "Pz", "Pt",
492 "PxAtDCA", "PyAtDCA", "PzAtDCA", "PtAtDCA",
493 "PxUncorrected", "PyUncorrected", "PzUncorrected", "PtUncorrected",
494 "XUncorrected", "YUncorrected", "ZUncorrected",
495 "XatDCA", "YatDCA", "DCA",
589e3f71 496 "Eta", "Rapidity", "Charge", "RAtAbsEnd",
b201705a 497 "IPVx", "IPVy", "IPVz"};
589e3f71 498 TString leavesInt[kNvarInt] = {"MatchTrig", "IsMuon", "IsGhost", "LoCircuit", "PassPhysicsSelection", "NVtxContrib", "NspdTracklets", "IsPileupVertex"};
b201705a 499 TString leavesChar[kNvarChar] = {"FiredTrigClass"};
b201705a 500 TString leavesUInt[kNvarUInt] = {"BunchCrossNum", "OrbitNum", "PeriodNum", "RunNum"};
d51d0b71 501 TString leavesFloatMC[kNvarFloatMC] = {"PxMC", "PyMC", "PzMC", "PtMC", "EtaMC", "RapidityMC", "VxMC", "VyMC", "VzMC"};
b201705a 502 TString leavesIntMC[kNvarIntMC] = {"Pdg", "MotherType"};
503
504 if ( ! fTreeSingleMu ) fTreeSingleMu = new TTree("fTreeSingleMu", "Single Mu");
505 if ( ! fTreeSingleMuMC ) fTreeSingleMuMC = new TTree("fTreeSingleMuMC", "Single Mu MC");
506
507 TTree* currTree[2] = {fTreeSingleMu, fTreeSingleMuMC};
b201705a 508
509 for(Int_t itree=0; itree<2; itree++){
589e3f71 510 TParameter<Int_t>* par1 = new TParameter<Int_t>("fillTreeScaleDown",fFillTreeScaleDown);
511 TParameter<Int_t>* par2 = new TParameter<Int_t>("keepAllEvents",fKeepAll);
512 currTree[itree]->GetUserInfo()->Add(par1);
513 currTree[itree]->GetUserInfo()->Add(par2);
b201705a 514 for(Int_t ivar=0; ivar<kNvarFloat; ivar++){
515 currTree[itree]->Branch(leavesFloat[ivar].Data(), &fVarFloat[ivar], Form("%s/F", leavesFloat[ivar].Data()));
516 }
517 for(Int_t ivar=0; ivar<kNvarInt; ivar++){
518 currTree[itree]->Branch(leavesInt[ivar].Data(), &fVarInt[ivar], Form("%s/I", leavesInt[ivar].Data()));
519 }
520 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
b201705a 521 TString addString = leavesChar[ivar] + "/C";
522 currTree[itree]->Branch(leavesChar[ivar].Data(), fVarChar[ivar], addString.Data());
523 }
524 for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
525 currTree[itree]->Branch(leavesUInt[ivar].Data(), &fVarUInt[ivar], Form("%s/i", leavesUInt[ivar].Data()));
526 }
527 if ( itree==1 ){
528 for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
529 currTree[itree]->Branch(leavesFloatMC[ivar].Data(), &fVarFloatMC[ivar], Form("%s/F", leavesFloatMC[ivar].Data()));
530 }
531 for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
532 currTree[itree]->Branch(leavesIntMC[ivar].Data(), &fVarIntMC[ivar], Form("%s/I", leavesIntMC[ivar].Data()));
533 }
534 }
535 } // loop on trees
536 } // fillNTuple
55065f3f 537
538 PostData(1, fCFManager->GetParticleContainer());
539 PostData(2, fHistoList);
540 PostData(3, fHistoListQA);
541 if ( fMCEvent ) PostData(4, fHistoListMC);
aad6618e 542}
543
544//________________________________________________________________________
9728bcfd 545void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/)
aad6618e 546{
547 //
548 /// Main loop
549 /// Called for each event
550 //
551
b201705a 552 AliESDEvent* esdEvent = 0x0;
9728bcfd 553 AliAODEvent* aodEvent = 0x0;
9728bcfd 554
b201705a 555 esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
556 if ( ! esdEvent ){
b201705a 557 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
558 }
aad6618e 559
b201705a 560 if ( ! aodEvent && ! esdEvent ) {
561 AliError ("AOD or ESD event not found. Nothing done!");
aad6618e 562 return;
563 }
564
589e3f71 565 if ( ! fMCEvent && InputEvent()->GetEventType() != 7 ) return; // Run only on physics events!
566
567 TTree* currTree = ( fMCEvent ) ? fTreeSingleMuMC : fTreeSingleMu;
568
569 Bool_t fillCurrentEventTree = ( fFillTreeScaleDown == 0 ) ? kFALSE : ( Entry() % fFillTreeScaleDown == 0 );
570
571 Reset(kFALSE);
572
573 // Global event info
574 TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses();
575 if ( fMCEvent ) {
576 // Add the MB name (which is not there in simulation)
577 // CAVEAT: to be checked if we perform beam-gas simulations
55065f3f 578 if ( ! firedTrigClasses.Contains("CINT") ) {
579 TString trigName = ((TObjString*)fTriggerClasses->At(0))->GetString();
580 firedTrigClasses.Prepend(Form("%s ", trigName.Data()));
581 }
589e3f71 582 }
55065f3f 583
589e3f71 584 fVarFloat[kVarIPVz] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetZ() : aodEvent->GetPrimaryVertex()->GetZ();
585 fVarInt[kVarNVtxContrib] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors();
586 fVarInt[kVarNspdTracklets] = ( esdEvent ) ? esdEvent->GetMultiplicity()->GetNumberOfTracklets() : aodEvent->GetTracklets()->GetNumberOfTracklets();
55065f3f 587 fVarInt[kVarIsPileup] = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8); // REMEMBER TO CHECK
589e3f71 588
589 fVarUInt[kVarRunNumber] = ( esdEvent ) ? esdEvent->GetRunNumber() : aodEvent->GetRunNumber();
590
591 Int_t isGoodVtxBin = ( fVarInt[kVarNVtxContrib] >= fkNvtxContribCut );
592 if ( fVarInt[kVarIsPileup] > 0 )
593 isGoodVtxBin = 2;
594
595 if ( fillCurrentEventTree ){
596 strcpy(fVarChar[kVarTrigMask], firedTrigClasses.Data());
b201705a 597
589e3f71 598 fVarInt[kVarPassPhysicsSelection] = ( esdEvent ) ? ((AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected() : ((AliAODInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
b201705a 599
600 // Small workaround: in MC the bunch ID are not properly set and the timestamp is in seconds
601 // So fill bunchCrossing with the read timestamp
602 // fill the orbit and period number with a timestamp created while reading the run
603 TTimeStamp ts;
589e3f71 604 fVarUInt[kVarBunchCrossNumber] = ( fMCEvent ) ? (UInt_t)Entry() : esdEvent->GetBunchCrossNumber();
605 fVarUInt[kVarOrbitNumber] = ( fMCEvent ) ? ts.GetNanoSec() : esdEvent->GetOrbitNumber();
606 fVarUInt[kVarPeriodNumber] = ( fMCEvent ) ? ts.GetTime() : esdEvent->GetPeriodNumber();
9728bcfd 607
589e3f71 608 fVarFloat[kVarIPVx] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetX() : aodEvent->GetPrimaryVertex()->GetX();
609 fVarFloat[kVarIPVy] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetY() : aodEvent->GetPrimaryVertex()->GetY();
610 }
aad6618e 611
612 // Object declaration
589e3f71 613 AliMCParticle* mcPart = 0x0;
614 AliVParticle* track = 0x0;
aad6618e 615
589e3f71 616 Double_t containerInput[kNvars];
617 Int_t histoIndex = -1;
b201705a 618
589e3f71 619 fCFManager->SetRecEventInfo(InputEvent());
b201705a 620
589e3f71 621 // Pure Monte Carlo part
622 if ( fMCEvent ) {
623 fCFManager->SetMCEventInfo (fMCEvent);
624 Int_t nMCtracks = fMCEvent->GetNumberOfTracks();
625 if ( nMCtracks > 0 ) {
626 containerInput[kHvarVz] = fMCEvent->Stack()->Particle(0)->Vz();
627 containerInput[kHvarIsGoodVtx] = 1;
628 histoIndex = GetHistoIndex(kHistoCheckVzMC);
629 ((TH2F*)fHistoListMC->At(histoIndex))->Fill(fVarFloat[kVarIPVz], fMCEvent->Stack()->Particle(0)->Vz());
630
631 }
632
633 for (Int_t ipart=0; ipart<nMCtracks; ipart++) {
634 mcPart = (AliMCParticle*)fMCEvent->GetTrack(ipart);
635
636 //check the MC-level cuts
637 if ( ! fCFManager->CheckParticleCuts(kStepGeneratedMC,mcPart) )
638 continue;
639
640 containerInput[kHvarPt] = mcPart->Pt();
55065f3f 641 containerInput[kHvarY] = mcPart->Y();
589e3f71 642 containerInput[kHvarPhi] = mcPart->Phi();
643 containerInput[kHvarDCA] = TMath::Sqrt(mcPart->Xv()*mcPart->Xv() +
644 mcPart->Yv()*mcPart->Yv());
645 containerInput[kHvarThetaZones] = GetBinThetaAbsEnd(TMath::Pi()-mcPart->Theta(),kTRUE);
646 containerInput[kHvarCharge] = mcPart->Charge()/3.;
647 containerInput[kHvarMatchTrig] = 1.;
589e3f71 648 containerInput[kHvarMotherType] = (Double_t)RecoTrackMother(mcPart);
649
55065f3f 650 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
651 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
652 if ( ! firedTrigClasses.Contains(trigName.Data()) )
653 continue;
654 containerInput[kHvarTrigClass] = (Double_t)(itrig+1);
655 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGeneratedMC);
656 if ( fCFManager->CheckParticleCuts(kStepAcceptanceMC,mcPart) )
657 fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptanceMC);
658 } // loop on trigger classes
589e3f71 659 } // loop on MC particles
660 } // is MC
661
662
663 Int_t trackLabel = -1;
b201705a 664 Bool_t isGhost = kFALSE;
665 Int_t nGhosts = 0, nMuons = 0;
666
589e3f71 667 Int_t nTracks = ( esdEvent ) ? esdEvent->GetNumberOfMuonTracks() : aodEvent->GetNTracks();
9728bcfd 668
669 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
aad6618e 670
9728bcfd 671 // Get variables
b201705a 672 if ( esdEvent ){
589e3f71 673 // ESD only info
9728bcfd 674
b201705a 675 track = esdEvent->GetMuonTrack(itrack);
676 isGhost = ( ((AliESDMuonTrack*)track)->ContainTriggerData() && ! ((AliESDMuonTrack*)track)->ContainTrackerData() );
677
589e3f71 678 if ( fillCurrentEventTree ) {
679 if ( itrack > 0 ) Reset(kTRUE);
680 fVarFloat[kVarPxUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaXUncorrected()) : ((AliESDMuonTrack*)track)->PxUncorrected();
681 fVarFloat[kVarPyUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaYUncorrected()) : ((AliESDMuonTrack*)track)->PyUncorrected();
682 fVarFloat[kVarPzUncorrected] = ( isGhost ) ? -1 : ((AliESDMuonTrack*)track)->PzUncorrected();
683
684 fVarFloat[kVarPtUncorrected] =
685 TMath::Sqrt(fVarFloat[kVarPxUncorrected] * fVarFloat[kVarPxUncorrected] +
686 fVarFloat[kVarPyUncorrected] * fVarFloat[kVarPyUncorrected]);
b201705a 687
589e3f71 688 fVarFloat[kVarXUncorrected] = ((AliESDMuonTrack*)track)->GetNonBendingCoorUncorrected();
689 fVarFloat[kVarYUncorrected] = ((AliESDMuonTrack*)track)->GetBendingCoorUncorrected();
690 fVarFloat[kVarZUncorrected] = ((AliESDMuonTrack*)track)->GetZUncorrected();
b201705a 691
589e3f71 692 fVarInt[kVarLocalCircuit] = ((AliESDMuonTrack*)track)->LoCircuit();
693 }
694
695 if ( isGhost ) {
696 // If is ghost fill only a partial information
697 nGhosts++;
698 if ( fillCurrentEventTree ) {
699 fVarInt[kVarIsMuon] = 0;
700 fVarInt[kVarIsGhost] = nGhosts;
701 fVarInt[kVarMatchTrig] = ((AliESDMuonTrack*)track)->GetMatchTrigger();
702 currTree->Fill();
703 }
b201705a 704 continue;
705 }
706 }
707 else {
708 track = aodEvent->GetTrack(itrack);
709 if ( ! ((AliAODTrack*)track)->IsMuonTrack() )
710 continue;
b201705a 711 }
9728bcfd 712
b201705a 713 // Information for tracks in tracker
714 nMuons++;
b201705a 715
716 fVarFloat[kVarPt] = track->Pt();
55065f3f 717 fVarFloat[kVarRapidity] = track->Y();
b201705a 718 fVarFloat[kVarXatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetNonBendingCoorAtDCA() : ((AliAODTrack*)track)->XAtDCA();
719 fVarFloat[kVarYatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetBendingCoorAtDCA() : ((AliAODTrack*)track)->YAtDCA();
589e3f71 720 fVarFloat[kVarDCA] =
721 TMath::Sqrt( fVarFloat[kVarXatDCA] * fVarFloat[kVarXatDCA] +
722 fVarFloat[kVarYatDCA] * fVarFloat[kVarYatDCA] );
723 fVarFloat[kVarCharge] = (Float_t)track->Charge();
724 fVarFloat[kVarRAtAbsEnd] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd() : ((AliAODTrack*)track)->GetRAtAbsorberEnd();
b201705a 725 fVarInt[kVarMatchTrig] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetMatchTrigger() : ((AliAODTrack*)track)->GetMatchTrigger();
726
589e3f71 727 if ( fillCurrentEventTree ){
728 fVarInt[kVarIsMuon] = nMuons;
729 fVarInt[kVarIsGhost] = 0;
55065f3f 730 fVarFloat[kVarEta] = track->Eta();
b201705a 731 fVarFloat[kVarPx] = track->Px();
732 fVarFloat[kVarPy] = track->Py();
733 fVarFloat[kVarPz] = track->Pz();
589e3f71 734 fVarFloat[kVarPxAtDCA] = ( esdEvent ) ? ((AliESDMuonTrack*)track)->PxAtDCA() : ((AliAODTrack*)track)->PxAtDCA();
735 fVarFloat[kVarPyAtDCA] = ( esdEvent ) ? ((AliESDMuonTrack*)track)->PyAtDCA() : ((AliAODTrack*)track)->PyAtDCA();
736 fVarFloat[kVarPzAtDCA] = ( esdEvent ) ? ((AliESDMuonTrack*)track)->PzAtDCA() : ((AliAODTrack*)track)->PzAtDCA();
b201705a 737 fVarFloat[kVarPtAtDCA] = TMath::Sqrt(fVarFloat[kVarPxAtDCA]*fVarFloat[kVarPxAtDCA] + fVarFloat[kVarPyAtDCA]*fVarFloat[kVarPyAtDCA]);
b201705a 738 }
9728bcfd 739
589e3f71 740 fVarIntMC[kVarMotherType] = kUnknownPart;
741
9728bcfd 742 // Monte Carlo part
589e3f71 743 if ( fMCEvent ) {
9728bcfd 744
589e3f71 745 trackLabel = track->GetLabel();
9728bcfd 746
589e3f71 747 AliMCParticle* matchedMCTrack = 0x0;
9728bcfd 748
589e3f71 749 Int_t nMCtracks = fMCEvent->GetNumberOfTracks();
750 for(Int_t imctrack=0; imctrack<nMCtracks; imctrack++){
751 mcPart = (AliMCParticle*)fMCEvent->GetTrack(imctrack);
752 if ( trackLabel == mcPart->GetLabel() ) {
753 matchedMCTrack = mcPart;
754 break;
755 }
756 } // loop on MC tracks
757
758 if ( matchedMCTrack ) {
759 fVarIntMC[kVarMotherType] = RecoTrackMother(matchedMCTrack);
760 fVarFloatMC[kVarPtMC] = matchedMCTrack->Pt();
761 FillTriggerHistos(kHistoPtResolutionMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt] - fVarFloatMC[kVarPtMC]);
762 if ( fillCurrentEventTree ){
763 fVarFloatMC[kVarPxMC] = matchedMCTrack->Px();
764 fVarFloatMC[kVarPyMC] = matchedMCTrack->Py();
765 fVarFloatMC[kVarPzMC] = matchedMCTrack->Pz();
766 fVarFloatMC[kVarEtaMC] = matchedMCTrack->Eta();
767 fVarFloatMC[kVarRapidityMC] = matchedMCTrack->Y();
768 fVarFloatMC[kVarVxMC] = matchedMCTrack->Xv();
769 fVarFloatMC[kVarVyMC] = matchedMCTrack->Yv();
770 fVarFloatMC[kVarVzMC] = matchedMCTrack->Zv();
771 fVarIntMC[kVarPdg] = matchedMCTrack->PdgCode();
772 }
b201705a 773 }
589e3f71 774 AliDebug(1, Form("Found mother %i", fVarIntMC[kVarMotherType]));
775 } // if use MC
776
777 containerInput[kHvarPt] = fVarFloat[kVarPt];
55065f3f 778 containerInput[kHvarY] = fVarFloat[kVarRapidity];
589e3f71 779 containerInput[kHvarPhi] = track->Phi();
780 containerInput[kHvarDCA] = fVarFloat[kVarDCA];
781 containerInput[kHvarVz] = fVarFloat[kVarIPVz];
782 containerInput[kHvarThetaZones] = GetBinThetaAbsEnd(fVarFloat[kVarRAtAbsEnd]);
783 containerInput[kHvarCharge] = fVarFloat[kVarCharge];
784 containerInput[kHvarMatchTrig] = (Double_t)fVarInt[kVarMatchTrig];
589e3f71 785 containerInput[kHvarIsGoodVtx] = (Double_t)isGoodVtxBin;
786 containerInput[kHvarMotherType] = (Double_t)fVarIntMC[kVarMotherType];
787
55065f3f 788 histoIndex = GetHistoIndex(kHistoNmuonsPerRun);
789 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
790 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
791 if ( ! firedTrigClasses.Contains(trigName.Data()) )
792 continue;
793 containerInput[kHvarTrigClass] = (Double_t)(itrig+1);
794 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed);
795 if ( fCFManager->CheckParticleCuts(kStepAcceptance,track) )
796 fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptance);
797 ((TH2F*)fHistoList->At(histoIndex))->Fill(Form("%u",fVarUInt[kVarRunNumber]), containerInput[kHvarTrigClass], 1.);
798 } // loop on trigger classes
589e3f71 799
800 if ( fillCurrentEventTree ) currTree->Fill();
9728bcfd 801 } // loop on tracks
b201705a 802
589e3f71 803 if ( fillCurrentEventTree && fKeepAll && ( ( nMuons + nGhosts ) == 0 ) ) {
b201705a 804 // Fill event also if there is not muon (when explicitely required)
589e3f71 805 currTree->Fill();
b201705a 806 }
807
55065f3f 808 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
809 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
810 if ( ! firedTrigClasses.Contains(trigName.Data()) )
811 continue;
812 Double_t trigClassBin = (Double_t)(itrig+1);
813 histoIndex = GetHistoIndex(kHistoNeventsPerTrig);
814 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 1., 1.); // All events
815 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 4., (Double_t)(fVarInt[kVarIsPileup]+1)); // All vertices
816 if ( isGoodVtxBin == 1 )
817 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 2., 1.); // Good vtx events
818 else if ( isGoodVtxBin == 2 ) {
819 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 3., 1.); // Pileup events
820 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 5., 1.); // Pileup vertices
821 }
589e3f71 822
55065f3f 823 histoIndex = GetHistoIndex(kHistoMuonMultiplicity);
824 ((TH1F*)fHistoList->At(histoIndex))->Fill(nMuons, trigClassBin);
589e3f71 825
55065f3f 826 if ( isGoodVtxBin == 1 && itrig == 0 ) {
827 histoIndex = GetHistoIndex(kHistoEventVz);
828 ((TH1F*)fHistoList->At(histoIndex))->Fill(fVarFloat[kVarIPVz]);
829
830 TString histoName = Form("hVtx%u",fVarUInt[kVarRunNumber]);
831 TH1F* histoVtxCurrRun = (TH1F*)fAuxObjects->GetValue(histoName.Data());
832 histoVtxCurrRun->Fill(fVarFloat[kVarIPVz]);
833 if ( nMuons + nGhosts > 0 )
834 histoVtxCurrRun->SetNormFactor(histoVtxCurrRun->GetNormFactor()+1.);
835 }
836
837 histoIndex = GetHistoIndex(kHistoNeventsPerRun);
838 ((TH2F*)fHistoList->At(histoIndex))->Fill(Form("%u",fVarUInt[kVarRunNumber]), trigClassBin, 1.);
839 } // loop on trigger classes
589e3f71 840
589e3f71 841
842 // Post final data. It will be written to a file with option "RECREATE"
843 PostData(1, fCFManager->GetParticleContainer());
844 PostData(2, fHistoList);
55065f3f 845 PostData(3, fHistoListQA);
846 if ( fMCEvent ) PostData(4, fHistoListMC);
589e3f71 847 if ( fFillTreeScaleDown > 0 )
848 PostData(5, currTree);
aad6618e 849}
850
851//________________________________________________________________________
852void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
853 //
854 /// Draw some histogram at the end.
855 //
589e3f71 856
857 AliCFContainer* container = dynamic_cast<AliCFContainer*> (GetOutputData(1));
858 if ( ! container ) {
859 AliError("Cannot find container in file");
860 return;
861 }
862
863 //Int_t histoIndex = -1;
864
865 if ( ! gROOT->IsBatch() ) {
93fd3322 866 TCanvas *c1_SingleMu = new TCanvas("c1_SingleMu","Vz vs Pt",10,10,310,310);
867 c1_SingleMu->SetFillColor(10); c1_SingleMu->SetHighLightColor(10);
868 c1_SingleMu->SetLeftMargin(0.15); c1_SingleMu->SetBottomMargin(0.15);
589e3f71 869 container->Project(kHvarPt,kHvarVz,kStepReconstructed)->Draw("COLZ");
aad6618e 870 }
589e3f71 871
aad6618e 872}
873
874//________________________________________________________________________
9728bcfd 875 void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
876 Float_t var1, Float_t var2, Float_t var3)
aad6618e 877{
878 //
9728bcfd 879 /// Fill all histograms passing the trigger cuts
aad6618e 880 //
662e37fe 881
9728bcfd 882 Int_t nTrigs = TMath::Max(1, matchTrig);
883 TArrayI trigToFill(nTrigs);
884 trigToFill[0] = matchTrig;
885 for(Int_t itrig = 1; itrig < matchTrig; itrig++){
886 trigToFill[itrig] = itrig;
887 }
662e37fe 888
9728bcfd 889 TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;
890
891 TString className;
892 for(Int_t itrig = 0; itrig < nTrigs; itrig++){
893 Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
894 className = histoList->At(currIndex)->ClassName();
895 AliDebug(3, Form("matchTrig %i Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
896 if (className.Contains("1"))
897 ((TH1F*)histoList->At(currIndex))->Fill(var1);
898 else if (className.Contains("2"))
899 ((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
900 else if (className.Contains("3"))
901 ((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
902 else
903 AliWarning("Histogram not filled");
904 }
aad6618e 905}
906
aad6618e 907//________________________________________________________________________
589e3f71 908Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcParticle)
aad6618e 909{
910 //
9728bcfd 911 /// Find track mother from kinematics
aad6618e 912 //
913
589e3f71 914 Int_t recoPdg = mcParticle->PdgCode();
aad6618e 915
9728bcfd 916 // Track is not a muon
589e3f71 917 if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
aad6618e 918
589e3f71 919 Int_t imother = mcParticle->GetMother();
aad6618e 920
589e3f71 921 Bool_t isFirstMotherHF = kFALSE;
922 Int_t step = 0;
aad6618e 923
589e3f71 924 while ( imother >= 0 ) {
925 TParticle* part = ((AliMCParticle*)fMCEvent->GetTrack(imother))->Particle();
662e37fe 926
589e3f71 927 if ( part->GetUniqueID() == kPHadronic )
928 return kSecondaryMu;
662e37fe 929
589e3f71 930 Int_t absPdg = TMath::Abs(part->GetPdgCode());
9728bcfd 931
589e3f71 932 step++;
933 if ( step == 1 ) // only for first mother
934 isFirstMotherHF = ( ( absPdg >= 400 && absPdg < 600 ) ||
935 ( absPdg >= 4000 && absPdg < 6000 ) );
936
937 if ( absPdg < 4 )
938 return kPrimaryMu;
9728bcfd 939
589e3f71 940 if ( isFirstMotherHF) {
941 if ( absPdg == 4 )
942 return kCharmMu;
943 else if ( absPdg == 5 )
944 return kBeautyMu;
9728bcfd 945 }
946
589e3f71 947 imother = part->GetFirstMother();
9728bcfd 948 }
949
589e3f71 950 return kPrimaryMu;
951}
662e37fe 952
589e3f71 953//________________________________________________________________________
954Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
955{
956 //
957 /// Get histogram index in the list
958 //
959
960 if ( srcIndex < 0 ) {
961 return histoTypeIndex;
9728bcfd 962 }
589e3f71 963
964 return
965 kNsummaryHistosMC +
966 kNtrackSources * kNtrigCuts * histoTypeIndex +
967 kNtrackSources * trigIndex +
968 srcIndex;
9728bcfd 969}
662e37fe 970
9728bcfd 971//________________________________________________________________________
589e3f71 972Float_t AliAnalysisTaskSingleMu::GetBinThetaAbsEnd(Float_t RAtAbsEnd, Bool_t isTheta)
973{
974 //
975 /// Get bin of theta at absorber end region
976 //
977 Float_t thetaDeg = ( isTheta ) ? RAtAbsEnd : TMath::ATan( RAtAbsEnd / 505. );
978 thetaDeg *= TMath::RadToDeg();
979 if ( thetaDeg < 2. )
980 return 0.;
981 else if ( thetaDeg < 3. )
982 return 1.;
983 else if ( thetaDeg < 9. )
984 return 2.;
985
986 return 3.;
987}
988
b201705a 989
990//________________________________________________________________________
991void AliAnalysisTaskSingleMu::Reset(Bool_t keepGlobal)
992{
993 //
994 /// Reset variables
995 //
996 Int_t lastVarFloat = ( keepGlobal ) ? kVarIPVx : kNvarFloat;
997 for(Int_t ivar=0; ivar<lastVarFloat; ivar++){
998 fVarFloat[ivar] = 0.;
999 }
1000 Int_t lastVarInt = ( keepGlobal ) ? kVarPassPhysicsSelection : kNvarInt;
1001 for(Int_t ivar=0; ivar<lastVarInt; ivar++){
589e3f71 1002 fVarInt[ivar] = 0;
b201705a 1003 }
589e3f71 1004 fVarInt[kVarMatchTrig] = -1;
b201705a 1005
1006 if ( ! keepGlobal ){
1007 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
1008 //sprintf(fVarChar[ivar],"%253s","");
1009 sprintf(fVarChar[ivar]," ");
1010 }
1011 for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
1012 fVarUInt[ivar] = 0;
1013 }
1014 }
589e3f71 1015 if ( fMCEvent ){
b201705a 1016 for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
1017 fVarFloatMC[ivar] = 0.;
1018 }
1019 for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
1020 fVarIntMC[ivar] = -1;
1021 }
1022 }
1023}
589e3f71 1024
1025//________________________________________________________________________
55065f3f 1026void AliAnalysisTaskSingleMu::SetTriggerClasses(TString triggerClasses)
589e3f71 1027{
55065f3f 1028 /// Set trigger classes
1029 if ( fTriggerClasses )
1030 delete fTriggerClasses;
1031
1032 fTriggerClasses = triggerClasses.Tokenize(" ");
1033 fTriggerClasses->SetOwner();
589e3f71 1034
55065f3f 1035}
589e3f71 1036
1037//________________________________________________________________________
1038void AliAnalysisTaskSingleMu::SetAxisLabel(TAxis* axis)
1039{
1040 //
1041 /// Set the labels of trigger class axis
1042 //
55065f3f 1043 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
1044 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
1045 axis->SetBinLabel(itrig+1,trigName.Data());
1046 }
1047 axis->SetTitle("Fired class");
589e3f71 1048}