]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/muon/AliAnalysisTaskSingleMu.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New function to get the list of keyWords of a given rubric (Philippe P.)
[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"
aad6618e 45
46// STEER includes
47#include "AliLog.h"
48
49#include "AliAODEvent.h"
50#include "AliAODTrack.h"
51#include "AliAODVertex.h"
aad6618e 52
9728bcfd 53#include "AliMCEvent.h"
54#include "AliMCParticle.h"
55#include "AliStack.h"
56
b201705a 57#include "AliESDInputHandler.h"
58#include "AliESDEvent.h"
59#include "AliESDMuonTrack.h"
60#include "AliAnalysisManager.h"
9728bcfd 61
62#include "AliAnalysisTaskSE.h"
aad6618e 63#include "AliAnalysisTaskSingleMu.h"
64
662e37fe 65/// \cond CLASSIMP
66ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
67/// \endcond
aad6618e 68
9728bcfd 69
aad6618e 70//________________________________________________________________________
b201705a 71AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, Bool_t fillTree, Bool_t keepAll) :
9728bcfd 72 AliAnalysisTaskSE(name),
73 fUseMC(0),
b201705a 74 fFillTree(fillTree),
75 fKeepAll(keepAll),
9728bcfd 76 fHistoList(0),
b201705a 77 fHistoListMC(0),
78 fTreeSingleMu(0),
79 fTreeSingleMuMC(0),
80 fVarFloat(0),
81 fVarInt(0),
82 fVarChar(0),
83 fVarUInt(0),
84 fVarFloatMC(0),
85 fVarIntMC(0)
aad6618e 86{
87 //
88 /// Constructor.
89 //
b201705a 90 if ( ! fFillTree )
91 fKeepAll = kFALSE;
92
9728bcfd 93 DefineOutput(1, TList::Class());
94 DefineOutput(2, TList::Class());
b201705a 95
96 if ( fFillTree )
97 DefineOutput(3, TTree::Class());
aad6618e 98}
99
9728bcfd 100
101//________________________________________________________________________
102AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
aad6618e 103{
104 //
9728bcfd 105 /// Destructor
aad6618e 106 //
b201705a 107 delete fHistoList;
108 delete fHistoListMC;
109 delete fTreeSingleMu;
110 delete fTreeSingleMuMC;
111 delete fVarFloat;
112 delete fVarInt;
113 delete [] fVarChar;
114 delete fVarUInt;
115 delete fVarFloatMC;
116 delete fVarIntMC;
9728bcfd 117}
aad6618e 118
9728bcfd 119
120//___________________________________________________________________________
b201705a 121void AliAnalysisTaskSingleMu::NotifyRun()
9728bcfd 122{
123 //
124 /// Use the event handler information to correctly fill the analysis flags:
125 /// - check if Monte Carlo information is present
126 //
127
128 if ( MCEvent() ) {
129 fUseMC = kTRUE;
130 AliInfo("Monte Carlo information is present");
131 }
132 else {
133 AliInfo("No Monte Carlo information in run");
aad6618e 134 }
135}
136
9728bcfd 137
aad6618e 138//___________________________________________________________________________
9728bcfd 139void AliAnalysisTaskSingleMu::UserCreateOutputObjects()
aad6618e 140{
141 //
142 /// Create output histograms
143 //
9728bcfd 144 AliInfo(Form(" CreateOutputObjects of task %s\n", GetName()));
aad6618e 145
b201705a 146 if ( fFillTree ) OpenFile(1);
147
9728bcfd 148 // initialize histogram lists
149 if(!fHistoList) fHistoList = new TList();
150 if(!fHistoListMC) fHistoListMC = new TList();
b201705a 151
152 // Init variables
153 fVarFloat = new Float_t [kNvarFloat];
154 fVarInt = new Int_t [kNvarInt];
155 fVarChar = new Char_t *[kNvarChar];
156 fVarUInt = new UInt_t [kNvarUInt];
157 fVarFloatMC = new Float_t [kNvarFloatMC];
158 fVarIntMC = new Int_t [kNvarIntMC];
9728bcfd 159
160 Int_t nPtBins = 60;
161 Float_t ptMin = 0., ptMax = 30.;
162 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
163
164 Int_t nDcaBins = 100;
165 Float_t dcaMin = 0., dcaMax = 200.;
166 TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");
662e37fe 167
9728bcfd 168 Int_t nVzBins = 60;
169 Float_t vzMin = -30, vzMax = 30;
170 TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
171
172 Int_t nEtaBins = 25;
173 Float_t etaMin = -4.5, etaMax = -2.;
174 TString etaName("Eta"), etaTitle("#eta"), etaUnits("a.u.");
175
176 Int_t nRapidityBins = 25;
177 Float_t rapidityMin = -4.5, rapidityMax = -2.;
178 TString rapidityName("Rapidity"), rapidityTitle("rapidity"), rapidityUnits("a.u.");
179
180 TString trigName[kNtrigCuts];
181 trigName[kNoMatchTrig] = "NoMatch";
182 trigName[kAllPtTrig] = "AllPt";
183 trigName[kLowPtTrig] = "LowPt";
184 trigName[kHighPtTrig] = "HighPt";
185
186 TString srcName[kNtrackSources];
187 srcName[kCharmMu] = "Charm";
188 srcName[kBeautyMu] = "Beauty";
189 srcName[kPrimaryMu] = "Decay";
190 srcName[kSecondaryMu] = "Secondary";
191 srcName[kRecoHadron] = "Hadrons";
192 srcName[kUnknownPart] = "Unknown";
193
194 TH1F* histo1D = 0x0;
195 TH2F* histo2D = 0x0;
196 TString histoName, histoTitle;
197 Int_t histoIndex = 0;
198
199 // 1D histos
200 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
201 histoName = Form("%sDistrib%sTrig", ptName.Data(), trigName[itrig].Data());
202 histoTitle = Form("%s distribution. Trigger: %s", ptTitle.Data(), trigName[itrig].Data());
203 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nPtBins, ptMin, ptMax);
204 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
205 histoIndex = GetHistoIndex(kHistoPt, itrig);
206 fHistoList->AddAt(histo1D, histoIndex);
207 }
208
209 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
210 histoName = Form("%sDistrib%sTrig", dcaName.Data(), trigName[itrig].Data());
211 histoTitle = Form("%s distribution. Trigger: %s", dcaTitle.Data(), trigName[itrig].Data());
212 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nDcaBins, dcaMin, dcaMax);
213 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
214 histoIndex = GetHistoIndex(kHistoDCA, itrig);
215 fHistoList->AddAt(histo1D, histoIndex);
216 }
217
218 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
219 histoName = Form("%sDistrib%sTrig", vzName.Data(), trigName[itrig].Data());
220 histoTitle = Form("%s distribution. Trigger: %s", vzTitle.Data(), trigName[itrig].Data());
221 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nVzBins, vzMin, vzMax);
222 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
223 histoIndex = GetHistoIndex(kHistoVz, itrig);
224 fHistoList->AddAt(histo1D, histoIndex);
225 }
226
227 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
228 histoName = Form("%sDistrib%sTrig", etaName.Data(), trigName[itrig].Data());
229 histoTitle = Form("%s distribution. Trigger: %s", etaTitle.Data(), trigName[itrig].Data());
230 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nEtaBins, etaMin, etaMax);
231 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", etaTitle.Data(), etaUnits.Data()));
232 histoIndex = GetHistoIndex(kHistoEta, itrig);
233 fHistoList->AddAt(histo1D, histoIndex);
234 }
235
236 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
237 histoName = Form("%sDistrib%sTrig", rapidityName.Data(), trigName[itrig].Data());
238 histoTitle = Form("%s distribution. Trigger: %s", rapidityTitle.Data(), trigName[itrig].Data());
239 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nRapidityBins, rapidityMin, rapidityMax);
240 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
241 histoIndex = GetHistoIndex(kHistoRapidity, itrig);
242 fHistoList->AddAt(histo1D, histoIndex);
243 }
662e37fe 244
9728bcfd 245 // 2D histos
246 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
247 histoName = Form("%s%sDistrib%sTrig", dcaName.Data(), ptName.Data(), trigName[itrig].Data());
248 histoTitle = Form("%s vs. %s distribution. Trigger: %s", dcaTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
249 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
250 nPtBins, ptMin, ptMax,
251 nDcaBins, dcaMin, dcaMax);
252 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
253 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
254 histoIndex = GetHistoIndex(kHistoPtDCA, itrig);
255 fHistoList->AddAt(histo2D, histoIndex);
256 }
257
258 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
259 histoName = Form("%s%sDistrib%sTrig", vzName.Data(), ptName.Data(), trigName[itrig].Data());
260 histoTitle = Form("%s vs. %s distribution. Trigger: %s", vzTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
261 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
262 nPtBins, ptMin, ptMax,
263 nVzBins, vzMin, vzMax);
264 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
265 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
266 histoIndex = GetHistoIndex(kHistoPtVz, itrig);
267 fHistoList->AddAt(histo2D, histoIndex);
268 }
269
270 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
271 histoName = Form("%s%sDistrib%sTrig", rapidityName.Data(), ptName.Data(), trigName[itrig].Data());
272 histoTitle = Form("%s vs. %s distribution. Trigger: %s", rapidityTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
273 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
274 nPtBins, ptMin, ptMax,
275 nRapidityBins, rapidityMin, rapidityMax);
276 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
277 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
278 histoIndex = GetHistoIndex(kHistoPtRapidity, itrig);
279 fHistoList->AddAt(histo2D, histoIndex);
662e37fe 280 }
281
b201705a 282 // Summary histos
283 histoName = "histoMuonMultiplicity";
284 histoTitle = "Muon track multiplicity";
285 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 10, -0.5, 10-0.5);
286 //histo1D->GetXaxis()->SetBinLabel(1, "All events");
287 //histo1D->GetXaxis()->SetBinLabel(2, "Muon events");
288 histo1D->GetXaxis()->SetTitle("# of muons");
289 histoIndex = GetHistoIndex(kNhistoTypes, 0) + kHistoMuonMultiplicity;
290 fHistoList->AddAt(histo1D, histoIndex);
291
9728bcfd 292 // MC histos
293 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
294 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
295 histoName = Form("%sResolution%sTrig%s", ptName.Data(), trigName[itrig].Data(), srcName[isrc].Data());
296 histoTitle = Form("%s resolution. Trigger: %s (%s)", ptTitle.Data(), trigName[itrig].Data(), srcName[isrc].Data());
297 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 100, -5., 5.);
298 histo1D->GetXaxis()->SetTitle(Form("%s^{reco} - %s^{MC} (%s)", ptTitle.Data(), ptTitle.Data(), ptUnits.Data()));
b201705a 299 histoIndex = GetHistoIndex(kHistoPtResolutionMC, itrig, isrc);
9728bcfd 300 fHistoListMC->AddAt(histo1D, histoIndex);
301 }
302 }
303
304 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
305 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
306 histoName = Form("%s%sDistrib%sTrig%s", dcaName.Data(), ptName.Data(),
307 trigName[itrig].Data(), srcName[isrc].Data());
308 histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", dcaTitle.Data(), ptTitle.Data(),
309 trigName[itrig].Data(), srcName[isrc].Data());
310 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
311 nPtBins, ptMin, ptMax,
312 nDcaBins, dcaMin, dcaMax);
313 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
314 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
b201705a 315 histoIndex = GetHistoIndex(kHistoPtDCAMC, itrig, isrc);
9728bcfd 316 fHistoListMC->AddAt(histo2D, histoIndex);
317 }
318 }
319
320 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
321 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
322 histoName = Form("%s%sDistrib%sTrig%s", vzName.Data(), ptName.Data(),
323 trigName[itrig].Data(), srcName[isrc].Data());
324 histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", vzTitle.Data(), ptTitle.Data(),
325 trigName[itrig].Data(), srcName[isrc].Data());
326 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
327 nPtBins, ptMin, ptMax,
328 nVzBins, vzMin, vzMax);
329 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
330 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
b201705a 331 histoIndex = GetHistoIndex(kHistoPtVzMC, itrig, isrc);
9728bcfd 332 fHistoListMC->AddAt(histo2D, histoIndex);
333 }
aad6618e 334 }
b201705a 335
336 // Trees
337 if ( fFillTree ){
338 TString leavesFloat[kNvarFloat] = {"Px", "Py", "Pz", "Pt",
339 "PxAtDCA", "PyAtDCA", "PzAtDCA", "PtAtDCA",
340 "PxUncorrected", "PyUncorrected", "PzUncorrected", "PtUncorrected",
341 "XUncorrected", "YUncorrected", "ZUncorrected",
342 "XatDCA", "YatDCA", "DCA",
343 "Eta", "Rapidity", "Charge",
344 "IPVx", "IPVy", "IPVz"};
345 TString leavesInt[kNvarInt] = {"MatchTrig", "IsMuon", "IsGhost", "PassPhysicsSelection"};
346 TString leavesChar[kNvarChar] = {"FiredTrigClass"};
347 const Int_t charWidth[kNvarChar] = {255};
348 TString leavesUInt[kNvarUInt] = {"BunchCrossNum", "OrbitNum", "PeriodNum", "RunNum"};
d51d0b71 349 TString leavesFloatMC[kNvarFloatMC] = {"PxMC", "PyMC", "PzMC", "PtMC", "EtaMC", "RapidityMC", "VxMC", "VyMC", "VzMC"};
b201705a 350 TString leavesIntMC[kNvarIntMC] = {"Pdg", "MotherType"};
351
352 if ( ! fTreeSingleMu ) fTreeSingleMu = new TTree("fTreeSingleMu", "Single Mu");
353 if ( ! fTreeSingleMuMC ) fTreeSingleMuMC = new TTree("fTreeSingleMuMC", "Single Mu MC");
354
355 TTree* currTree[2] = {fTreeSingleMu, fTreeSingleMuMC};
356 //TList* currList[2] = {fHistoList, fHistoListMC};
357
358 for(Int_t itree=0; itree<2; itree++){
359 for(Int_t ivar=0; ivar<kNvarFloat; ivar++){
360 currTree[itree]->Branch(leavesFloat[ivar].Data(), &fVarFloat[ivar], Form("%s/F", leavesFloat[ivar].Data()));
361 }
362 for(Int_t ivar=0; ivar<kNvarInt; ivar++){
363 currTree[itree]->Branch(leavesInt[ivar].Data(), &fVarInt[ivar], Form("%s/I", leavesInt[ivar].Data()));
364 }
365 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
366 if ( itree == 0 ){
367 fVarChar[ivar] = new Char_t [charWidth[ivar]];
368 }
369 TString addString = leavesChar[ivar] + "/C";
370 currTree[itree]->Branch(leavesChar[ivar].Data(), fVarChar[ivar], addString.Data());
371 }
372 for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
373 currTree[itree]->Branch(leavesUInt[ivar].Data(), &fVarUInt[ivar], Form("%s/i", leavesUInt[ivar].Data()));
374 }
375 if ( itree==1 ){
376 for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
377 currTree[itree]->Branch(leavesFloatMC[ivar].Data(), &fVarFloatMC[ivar], Form("%s/F", leavesFloatMC[ivar].Data()));
378 }
379 for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
380 currTree[itree]->Branch(leavesIntMC[ivar].Data(), &fVarIntMC[ivar], Form("%s/I", leavesIntMC[ivar].Data()));
381 }
382 }
383 } // loop on trees
384 } // fillNTuple
aad6618e 385}
386
387//________________________________________________________________________
9728bcfd 388void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/)
aad6618e 389{
390 //
391 /// Main loop
392 /// Called for each event
393 //
394
b201705a 395 AliESDEvent* esdEvent = 0x0;
9728bcfd 396 AliAODEvent* aodEvent = 0x0;
397 AliMCEvent* mcEvent = 0x0;
398
b201705a 399 esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
400 if ( ! esdEvent ){
401 fFillTree = kFALSE;
402 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
403 }
404 //else
405 //aodEvent = AODEvent();
aad6618e 406
b201705a 407 if ( ! aodEvent && ! esdEvent ) {
408 AliError ("AOD or ESD event not found. Nothing done!");
aad6618e 409 return;
410 }
411
b201705a 412 if ( ! fUseMC && InputEvent()->GetEventType() !=7 ) return; // Run only on physics events!
413
414 if ( fFillTree ){
415 Reset(kFALSE);
416 fVarInt[kVarPassPhysicsSelection] = ((AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
417 strcpy(fVarChar[kVarTrigMask], esdEvent->GetFiredTriggerClasses().Data());
418
419 // Small workaround: in MC the bunch ID are not properly set and the timestamp is in seconds
420 // So fill bunchCrossing with the read timestamp
421 // fill the orbit and period number with a timestamp created while reading the run
422 TTimeStamp ts;
d51d0b71 423 fVarUInt[kVarBunchCrossNumber] = ( fUseMC ) ? (UInt_t)Entry() : esdEvent->GetBunchCrossNumber();
b201705a 424 fVarUInt[kVarOrbitNumber] = ( fUseMC ) ? ts.GetNanoSec() : esdEvent->GetOrbitNumber();
425 fVarUInt[kVarPeriodNumber] = ( fUseMC ) ? ts.GetTime() : esdEvent->GetPeriodNumber();
426 fVarUInt[kVarRunNumber] = esdEvent->GetRunNumber();
427
428 fVarFloat[kVarIPVx] = esdEvent->GetPrimaryVertex()->GetX();
429 fVarFloat[kVarIPVy] = esdEvent->GetPrimaryVertex()->GetY();
430 fVarFloat[kVarIPVz] = esdEvent->GetPrimaryVertex()->GetZ();
431 }
9728bcfd 432
433 if ( fUseMC ) mcEvent = MCEvent();
aad6618e 434
435 // Object declaration
9728bcfd 436 AliMCParticle* mcTrack = 0x0;
aad6618e 437
9728bcfd 438 Int_t trackLabel = -1;
b201705a 439
440 Int_t nTracks = ( esdEvent ) ? esdEvent->GetNumberOfMuonTracks() : aodEvent->GetNTracks();
441
442 Bool_t isGhost = kFALSE;
443 Int_t nGhosts = 0, nMuons = 0;
444
445 AliVParticle* track = 0x0;
9728bcfd 446
447 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
aad6618e 448
9728bcfd 449 // Get variables
b201705a 450 if ( esdEvent ){
451 if (itrack>0) Reset(kTRUE);
9728bcfd 452
b201705a 453 track = esdEvent->GetMuonTrack(itrack);
454 isGhost = ( ((AliESDMuonTrack*)track)->ContainTriggerData() && ! ((AliESDMuonTrack*)track)->ContainTrackerData() );
455
456 // ESD only info
457 fVarFloat[kVarPxUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaXUncorrected()) : ((AliESDMuonTrack*)track)->PxUncorrected();
458 fVarFloat[kVarPyUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaYUncorrected()) : ((AliESDMuonTrack*)track)->PyUncorrected();
459 fVarFloat[kVarPzUncorrected] = ( isGhost ) ? -1 : ((AliESDMuonTrack*)track)->PzUncorrected();
460 fVarFloat[kVarPtUncorrected] = TMath::Sqrt(
461 fVarFloat[kVarPxUncorrected] * fVarFloat[kVarPxUncorrected] +
462 fVarFloat[kVarPyUncorrected] * fVarFloat[kVarPyUncorrected]);
463
464 fVarFloat[kVarXUncorrected] = ((AliESDMuonTrack*)track)->GetNonBendingCoorUncorrected();
465 fVarFloat[kVarYUncorrected] = ((AliESDMuonTrack*)track)->GetBendingCoorUncorrected();
466 fVarFloat[kVarZUncorrected] = ((AliESDMuonTrack*)track)->GetZUncorrected();
467
468 fVarInt[kVarMatchTrig] = ((AliESDMuonTrack*)track)->GetMatchTrigger();
469
470 // If is ghost fill only a partial information
471 if ( isGhost ){
472 fVarInt[kVarIsMuon] = 0;
473 nGhosts++;
474 fVarInt[kVarIsGhost] = nGhosts;
475
476 if ( ! fUseMC ) fTreeSingleMu->Fill();
477 else fTreeSingleMuMC->Fill();
478
479 continue;
480 }
481 }
482 else {
483 track = aodEvent->GetTrack(itrack);
484 if ( ! ((AliAODTrack*)track)->IsMuonTrack() )
485 continue;
486
487 //Reset(kTRUE);
488 }
9728bcfd 489
b201705a 490 // Information for tracks in tracker
491 nMuons++;
492 fVarInt[kVarIsMuon] = nMuons;
493 fVarInt[kVarIsGhost] = 0;
494
495 fVarFloat[kVarPt] = track->Pt();
496 fVarFloat[kVarXatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetNonBendingCoorAtDCA() : ((AliAODTrack*)track)->XAtDCA();
497 fVarFloat[kVarYatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetBendingCoorAtDCA() : ((AliAODTrack*)track)->YAtDCA();
498 fVarFloat[kVarDCA] = TMath::Sqrt( fVarFloat[kVarXatDCA] * fVarFloat[kVarXatDCA] + fVarFloat[kVarYatDCA] * fVarFloat[kVarYatDCA] );
499 fVarFloat[kVarEta] = track->Eta();
500 fVarFloat[kVarRapidity] = track->Y();
501 trackLabel = track->GetLabel();
502 fVarInt[kVarMatchTrig] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetMatchTrigger() : ((AliAODTrack*)track)->GetMatchTrigger();
503
504 // Fill histograms
505 FillTriggerHistos(kHistoPt, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt]);
506 FillTriggerHistos(kHistoDCA, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarDCA]);
507 FillTriggerHistos(kHistoVz, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarIPVz]);
508 FillTriggerHistos(kHistoEta, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarEta]);
509 FillTriggerHistos(kHistoRapidity, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarRapidity]);
510
511 FillTriggerHistos(kHistoPtDCA, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt], fVarFloat[kVarDCA]);
512 FillTriggerHistos(kHistoPtVz, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt], fVarFloat[kVarIPVz]);
513 FillTriggerHistos(kHistoPtRapidity, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt], fVarFloat[kVarRapidity]);
514
515 if ( fFillTree ){
516 fVarFloat[kVarPx] = track->Px();
517 fVarFloat[kVarPy] = track->Py();
518 fVarFloat[kVarPz] = track->Pz();
519 fVarFloat[kVarPxAtDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->PxAtDCA() : ((AliAODTrack*)track)->PxAtDCA();
520 fVarFloat[kVarPyAtDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->PyAtDCA() : ((AliAODTrack*)track)->PyAtDCA();
521 fVarFloat[kVarPzAtDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->PzAtDCA() : ((AliAODTrack*)track)->PzAtDCA();
522 fVarFloat[kVarPtAtDCA] = TMath::Sqrt(fVarFloat[kVarPxAtDCA]*fVarFloat[kVarPxAtDCA] + fVarFloat[kVarPyAtDCA]*fVarFloat[kVarPyAtDCA]);
523
524 fVarFloat[kVarCharge] = track->Charge();
525
526 if ( ! fUseMC ) fTreeSingleMu->Fill();
527 }
9728bcfd 528
529 // Monte Carlo part
b201705a 530 if ( ! fUseMC ) continue;
9728bcfd 531
b201705a 532 fVarIntMC[kVarMotherType] = kUnknownPart;
9728bcfd 533
534 AliMCParticle* matchedMCTrack = 0x0;
535
536 Int_t nMCtracks = mcEvent->GetNumberOfTracks();
537 for(Int_t imctrack=0; imctrack<nMCtracks; imctrack++){
538 mcTrack = (AliMCParticle*)mcEvent->GetTrack(imctrack);
539 if ( trackLabel == mcTrack->GetLabel() ) {
540 matchedMCTrack = mcTrack;
541 break;
542 }
543 } // loop on MC tracks
544
545 if ( matchedMCTrack )
b201705a 546 fVarIntMC[kVarMotherType] = RecoTrackMother(matchedMCTrack, mcEvent);
9728bcfd 547
b201705a 548 AliDebug(1, Form("Found mother %i", fVarIntMC[kVarMotherType]));
9728bcfd 549
9728bcfd 550
b201705a 551 FillTriggerHistos(kHistoPtDCAMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt], fVarFloat[kVarDCA]);
552 FillTriggerHistos(kHistoPtVzMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt], fVarFloat[kVarIPVz]);
553
d51d0b71 554 if ( matchedMCTrack ) {
555 fVarFloatMC[kVarPtMC] = matchedMCTrack->Pt();
b201705a 556 FillTriggerHistos(kHistoPtResolutionMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt] - fVarFloatMC[kVarPtMC]);
557 if ( fFillTree ){
d51d0b71 558 fVarFloatMC[kVarPxMC] = matchedMCTrack->Px();
559 fVarFloatMC[kVarPyMC] = matchedMCTrack->Py();
560 fVarFloatMC[kVarPzMC] = matchedMCTrack->Pz();
561 fVarFloatMC[kVarEtaMC] = matchedMCTrack->Eta();
562 fVarFloatMC[kVarRapidityMC] = matchedMCTrack->Y();
563 fVarFloatMC[kVarVxMC] = matchedMCTrack->Xv();
564 fVarFloatMC[kVarVyMC] = matchedMCTrack->Yv();
565 fVarFloatMC[kVarVzMC] = matchedMCTrack->Zv();
566 fVarIntMC[kVarPdg] = matchedMCTrack->PdgCode();
b201705a 567 }
568 }
d51d0b71 569 if ( fFillTree ) fTreeSingleMuMC->Fill();
9728bcfd 570 } // loop on tracks
b201705a 571
572 if ( fKeepAll && ( ( fVarInt[kVarIsMuon] + fVarInt[kVarIsGhost] ) == 0 ) ) {
573 // Fill event also if there is not muon (when explicitely required)
574 if ( ! fUseMC ) fTreeSingleMu->Fill();
575 else fTreeSingleMuMC->Fill();
576 }
577
578 if( strstr(fVarChar[kVarTrigMask],"MUON") && fVarInt[kVarIsMuon]==0 )
579 printf("WARNING: Muon trigger does not match tracker!\n");
580
581 Int_t histoIndex = GetHistoIndex(kNhistoTypes, 0) + kHistoMuonMultiplicity;
582 ((TH1F*)fHistoList->At(histoIndex))->Fill(fVarInt[kVarIsMuon]);
583 //if ( fVarInt[kVarIsMuon]>0 ) ((TH1F*)fHistoList->At(histoIndex))->Fill(2);
9728bcfd 584
aad6618e 585 // Post final data. It will be written to a file with option "RECREATE"
9728bcfd 586 PostData(1, fHistoList);
587 if ( fUseMC ) PostData(2, fHistoListMC);
b201705a 588 if ( fFillTree ){
589 if ( fUseMC )
590 PostData(3, fTreeSingleMuMC);
591 else
592 PostData(3, fTreeSingleMu);
593 }
aad6618e 594}
595
596//________________________________________________________________________
597void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
598 //
599 /// Draw some histogram at the end.
600 //
601 if (!gROOT->IsBatch()) {
93fd3322 602 fHistoList = dynamic_cast<TList*> (GetOutputData(1));
603 TCanvas *c1_SingleMu = new TCanvas("c1_SingleMu","Vz vs Pt",10,10,310,310);
604 c1_SingleMu->SetFillColor(10); c1_SingleMu->SetHighLightColor(10);
605 c1_SingleMu->SetLeftMargin(0.15); c1_SingleMu->SetBottomMargin(0.15);
9728bcfd 606 Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
607 ((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
aad6618e 608 }
609}
610
611//________________________________________________________________________
9728bcfd 612 void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
613 Float_t var1, Float_t var2, Float_t var3)
aad6618e 614{
615 //
9728bcfd 616 /// Fill all histograms passing the trigger cuts
aad6618e 617 //
662e37fe 618
9728bcfd 619 Int_t nTrigs = TMath::Max(1, matchTrig);
620 TArrayI trigToFill(nTrigs);
621 trigToFill[0] = matchTrig;
622 for(Int_t itrig = 1; itrig < matchTrig; itrig++){
623 trigToFill[itrig] = itrig;
624 }
662e37fe 625
9728bcfd 626 TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;
627
628 TString className;
629 for(Int_t itrig = 0; itrig < nTrigs; itrig++){
630 Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
631 className = histoList->At(currIndex)->ClassName();
632 AliDebug(3, Form("matchTrig %i Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
633 if (className.Contains("1"))
634 ((TH1F*)histoList->At(currIndex))->Fill(var1);
635 else if (className.Contains("2"))
636 ((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
637 else if (className.Contains("3"))
638 ((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
639 else
640 AliWarning("Histogram not filled");
641 }
aad6618e 642}
643
aad6618e 644//________________________________________________________________________
9728bcfd 645Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
aad6618e 646{
647 //
9728bcfd 648 /// Find track mother from kinematics
aad6618e 649 //
650
9728bcfd 651 Int_t recoPdg = mcTrack->PdgCode();
662e37fe 652
9728bcfd 653 Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;
aad6618e 654
9728bcfd 655 // Track is not a muon
656 if ( ! isMuon ) return kRecoHadron;
aad6618e 657
9728bcfd 658 Int_t imother = mcTrack->GetMother();
aad6618e 659
9728bcfd 660 if ( imother<0 )
661 return kPrimaryMu; // Drell-Yan Muon
aad6618e 662
9728bcfd 663 Int_t igrandma = imother;
662e37fe 664
9728bcfd 665 AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
666 Int_t motherPdg = motherPart->PdgCode();
662e37fe 667
9728bcfd 668 // Track is an heavy flavor muon
669 Int_t absPdg = TMath::Abs(motherPdg);
670 if(absPdg/100==5 || absPdg/1000==5) {
671 return kBeautyMu;
672 }
673 if(absPdg/100==4 || absPdg/1000==4){
674 Int_t newMother = -1;
675 igrandma = imother;
676 AliInfo("\nFound candidate B->C->mu. History:\n");
677 mcTrack->Print();
678 printf("- %2i ", igrandma);
679 motherPart->Print();
680 Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
681 while ( absGrandMotherPdg > 10 ) {
682 igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
683 if( igrandma < 0 ) break;
684 printf("- %2i ", igrandma);
685 mcEvent->GetTrack(igrandma)->Print();
686 absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
687 }
688
689 if (absGrandMotherPdg==5) newMother = kBeautyMu; // Charm from beauty
690 else if (absGrandMotherPdg==4) newMother = kCharmMu;
691
692 if(newMother<0) {
693 AliWarning("Mother not correctly found! Set to charm!\n");
694 newMother = kCharmMu;
695 }
696
697 return newMother;
698 }
699
700 Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
662e37fe 701
9728bcfd 702 // Track is a bkg. muon
703 if (imother<nPrimaries) {
704 return kPrimaryMu;
705 }
706 else {
707 return kSecondaryMu;
708 }
709}
662e37fe 710
9728bcfd 711//________________________________________________________________________
712Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
713 {
b201705a 714 //
715 /// Get histogram index in the list
716 //
9728bcfd 717 if ( srcIndex < 0 ) return kNtrigCuts * histoTypeIndex + trigIndex;
718
719 return
720 kNtrackSources * kNtrigCuts * histoTypeIndex +
721 kNtrackSources * trigIndex +
722 srcIndex;
aad6618e 723}
b201705a 724
725//________________________________________________________________________
726void AliAnalysisTaskSingleMu::Reset(Bool_t keepGlobal)
727{
728 //
729 /// Reset variables
730 //
731 Int_t lastVarFloat = ( keepGlobal ) ? kVarIPVx : kNvarFloat;
732 for(Int_t ivar=0; ivar<lastVarFloat; ivar++){
733 fVarFloat[ivar] = 0.;
734 }
735 Int_t lastVarInt = ( keepGlobal ) ? kVarPassPhysicsSelection : kNvarInt;
736 for(Int_t ivar=0; ivar<lastVarInt; ivar++){
737 fVarInt[ivar] = -1;
738 }
739 fVarInt[kVarIsMuon] = 0;
740 fVarInt[kVarIsGhost] = 0;
741
742 if ( ! keepGlobal ){
743 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
744 //sprintf(fVarChar[ivar],"%253s","");
745 sprintf(fVarChar[ivar]," ");
746 }
747 for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
748 fVarUInt[ivar] = 0;
749 }
750 }
751 if ( fUseMC ){
752 for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
753 fVarFloatMC[ivar] = 0.;
754 }
755 for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
756 fVarIntMC[ivar] = -1;
757 }
758 }
759}