]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG3/muon/AliAnalysisTaskSingleMu.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-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
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
16//-----------------------------------------------------------------------------
17/// \class AliAnalysisTaskSingleMu
18/// Analysis task for single muons in the spectrometer.
19/// The output is a list of histograms.
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.
22/// If Monte Carlo information is present, some basics checks are performed.
23///
24/// \author Diego Stocco
25//-----------------------------------------------------------------------------
26
27//----------------------------------------------------------------------------
28// Implementation of the single muon analysis class
29// An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
30//----------------------------------------------------------------------------
31
32#define AliAnalysisTaskSingleMu_cxx
33
34// ROOT includes
35#include "TROOT.h"
36#include "TH1.h"
37#include "TH2.h"
38#include "TH3.h"
39#include "TAxis.h"
40#include "TList.h"
41#include "TCanvas.h"
42#include "TMath.h"
43#include "TTree.h"
44#include "TTimeStamp.h"
45
46// STEER includes
47#include "AliLog.h"
48
49#include "AliAODEvent.h"
50#include "AliAODTrack.h"
51#include "AliAODVertex.h"
52
53#include "AliMCEvent.h"
54#include "AliMCParticle.h"
55#include "AliStack.h"
56
57#include "AliESDInputHandler.h"
58#include "AliESDEvent.h"
59#include "AliESDMuonTrack.h"
60#include "AliAnalysisManager.h"
61
62#include "AliAnalysisTaskSE.h"
63#include "AliAnalysisTaskSingleMu.h"
64
65/// \cond CLASSIMP
66ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
67/// \endcond
68
69
70//________________________________________________________________________
71AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, Bool_t fillTree, Bool_t keepAll) :
72 AliAnalysisTaskSE(name),
73 fUseMC(0),
74 fFillTree(fillTree),
75 fKeepAll(keepAll),
76 fHistoList(0),
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)
86{
87 //
88 /// Constructor.
89 //
90 if ( ! fFillTree )
91 fKeepAll = kFALSE;
92
93 DefineOutput(1, TList::Class());
94 DefineOutput(2, TList::Class());
95
96 if ( fFillTree )
97 DefineOutput(3, TTree::Class());
98}
99
100
101//________________________________________________________________________
102AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
103{
104 //
105 /// Destructor
106 //
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;
117}
118
119
120//___________________________________________________________________________
121void AliAnalysisTaskSingleMu::NotifyRun()
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");
134 }
135}
136
137
138//___________________________________________________________________________
139void AliAnalysisTaskSingleMu::UserCreateOutputObjects()
140{
141 //
142 /// Create output histograms
143 //
144 AliInfo(Form(" CreateOutputObjects of task %s\n", GetName()));
145
146 if ( fFillTree ) OpenFile(1);
147
148 // initialize histogram lists
149 if(!fHistoList) fHistoList = new TList();
150 if(!fHistoListMC) fHistoListMC = new TList();
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];
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");
167
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 }
244
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);
280 }
281
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
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()));
299 histoIndex = GetHistoIndex(kHistoPtResolutionMC, itrig, isrc);
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()));
315 histoIndex = GetHistoIndex(kHistoPtDCAMC, itrig, isrc);
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()));
331 histoIndex = GetHistoIndex(kHistoPtVzMC, itrig, isrc);
332 fHistoListMC->AddAt(histo2D, histoIndex);
333 }
334 }
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"};
349 TString leavesFloatMC[kNvarFloatMC] = {"PxMC", "PyMC", "PzMC", "PtMC", "EtaMC", "RapidityMC", "VxMC", "VyMC", "VzMC"};
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
385}
386
387//________________________________________________________________________
388void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/)
389{
390 //
391 /// Main loop
392 /// Called for each event
393 //
394
395 AliESDEvent* esdEvent = 0x0;
396 AliAODEvent* aodEvent = 0x0;
397 AliMCEvent* mcEvent = 0x0;
398
399 esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
400 if ( ! esdEvent ){
401 fFillTree = kFALSE;
402 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
403 }
404 //else
405 //aodEvent = AODEvent();
406
407 if ( ! aodEvent && ! esdEvent ) {
408 AliError ("AOD or ESD event not found. Nothing done!");
409 return;
410 }
411
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;
423 fVarUInt[kVarBunchCrossNumber] = ( fUseMC ) ? (UInt_t)Entry() : esdEvent->GetBunchCrossNumber();
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 }
432
433 if ( fUseMC ) mcEvent = MCEvent();
434
435 // Object declaration
436 AliMCParticle* mcTrack = 0x0;
437
438 Int_t trackLabel = -1;
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;
446
447 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
448
449 // Get variables
450 if ( esdEvent ){
451 if (itrack>0) Reset(kTRUE);
452
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 }
489
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 }
528
529 // Monte Carlo part
530 if ( ! fUseMC ) continue;
531
532 fVarIntMC[kVarMotherType] = kUnknownPart;
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 )
546 fVarIntMC[kVarMotherType] = RecoTrackMother(matchedMCTrack, mcEvent);
547
548 AliDebug(1, Form("Found mother %i", fVarIntMC[kVarMotherType]));
549
550
551 FillTriggerHistos(kHistoPtDCAMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt], fVarFloat[kVarDCA]);
552 FillTriggerHistos(kHistoPtVzMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt], fVarFloat[kVarIPVz]);
553
554 if ( matchedMCTrack ) {
555 fVarFloatMC[kVarPtMC] = matchedMCTrack->Pt();
556 FillTriggerHistos(kHistoPtResolutionMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt] - fVarFloatMC[kVarPtMC]);
557 if ( fFillTree ){
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();
567 }
568 }
569 if ( fFillTree ) fTreeSingleMuMC->Fill();
570 } // loop on tracks
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);
584
585 // Post final data. It will be written to a file with option "RECREATE"
586 PostData(1, fHistoList);
587 if ( fUseMC ) PostData(2, fHistoListMC);
588 if ( fFillTree ){
589 if ( fUseMC )
590 PostData(3, fTreeSingleMuMC);
591 else
592 PostData(3, fTreeSingleMu);
593 }
594}
595
596//________________________________________________________________________
597void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
598 //
599 /// Draw some histogram at the end.
600 //
601 if (!gROOT->IsBatch()) {
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);
606 Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
607 ((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
608 }
609}
610
611//________________________________________________________________________
612 void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
613 Float_t var1, Float_t var2, Float_t var3)
614{
615 //
616 /// Fill all histograms passing the trigger cuts
617 //
618
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 }
625
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 }
642}
643
644//________________________________________________________________________
645Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
646{
647 //
648 /// Find track mother from kinematics
649 //
650
651 Int_t recoPdg = mcTrack->PdgCode();
652
653 Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;
654
655 // Track is not a muon
656 if ( ! isMuon ) return kRecoHadron;
657
658 Int_t imother = mcTrack->GetMother();
659
660 if ( imother<0 )
661 return kPrimaryMu; // Drell-Yan Muon
662
663 Int_t igrandma = imother;
664
665 AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
666 Int_t motherPdg = motherPart->PdgCode();
667
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();
701
702 // Track is a bkg. muon
703 if (imother<nPrimaries) {
704 return kPrimaryMu;
705 }
706 else {
707 return kSecondaryMu;
708 }
709}
710
711//________________________________________________________________________
712Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
713 {
714 //
715 /// Get histogram index in the list
716 //
717 if ( srcIndex < 0 ) return kNtrigCuts * histoTypeIndex + trigIndex;
718
719 return
720 kNtrackSources * kNtrigCuts * histoTypeIndex +
721 kNtrackSources * trigIndex +
722 srcIndex;
723}
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}