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