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