Casting TH1 to TH2
[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskSingleMu.cxx
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 #include "TMap.h"
46 #include "TObjString.h"
47 #include "TObjArray.h"
48 #include "TIterator.h"
49 #include "TParameter.h"
50 #include "TMCProcess.h"
51
52 // STEER includes
53 #include "AliLog.h"
54
55 #include "AliAODInputHandler.h"
56 #include "AliAODEvent.h"
57 #include "AliAODTrack.h"
58 #include "AliAODVertex.h"
59
60 #include "AliMCParticle.h"
61 #include "AliStack.h"
62
63 #include "AliESDInputHandler.h"
64 #include "AliESDEvent.h"
65 #include "AliESDMuonTrack.h"
66
67 #include "AliMultiplicity.h"
68
69 // ANALYSIS includes
70 #include "AliAnalysisManager.h"
71 #include "AliAnalysisTaskSE.h"
72 #include "AliAnalysisDataSlot.h"
73 #include "AliAnalysisDataContainer.h"
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
83 #include "AliAnalysisTaskSingleMu.h"
84
85 /// \cond CLASSIMP
86 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
87 /// \endcond
88
89
90 //________________________________________________________________________
91 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, Int_t fillTreeScaleDown, Bool_t keepAll) :
92   AliAnalysisTaskSE(name),
93   fFillTreeScaleDown(fillTreeScaleDown),
94   fKeepAll(keepAll),
95   fkNvtxContribCut(1),
96   fTriggerClasses(0x0),
97   fCFManager(0),
98   fHistoList(0),
99   fHistoListMC(0),
100   fHistoListQA(0),
101   fTreeSingleMu(0),
102   fVarFloat(0),
103   fVarInt(0),
104   fVarChar(0),
105   fVarUInt(0),
106   fVarFloatMC(0),
107   fVarIntMC(0),
108   fAuxObjects(new TMap())
109 {
110   //
111   /// Constructor.
112   //
113   if ( fFillTreeScaleDown <= 0 )
114     fKeepAll = kFALSE;
115
116   DefineOutput(1, AliCFContainer::Class());
117   DefineOutput(2, TList::Class());
118   DefineOutput(3, TList::Class());
119   DefineOutput(4, TList::Class());
120
121   if ( fFillTreeScaleDown > 0 )
122     DefineOutput(5, TTree::Class());
123
124   fAuxObjects->SetOwner();
125
126   SetTriggerClasses();
127 }
128
129
130 //________________________________________________________________________
131 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
132 {
133   //
134   /// Destructor
135   //
136
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   }
147   delete fVarFloat;
148   delete fVarInt;
149   delete [] fVarChar;
150   delete fVarUInt;
151   delete fVarFloatMC;
152   delete fVarIntMC;
153   delete fAuxObjects;
154 }
155
156
157 //___________________________________________________________________________
158 void AliAnalysisTaskSingleMu::NotifyRun()
159 {
160   //
161   /// Use the event handler information to correctly fill the analysis flags:
162   /// - check if Monte Carlo information is present
163   //  
164   
165   if ( fMCEvent ) {
166     AliInfo("Monte Carlo information is present");
167   }
168   else {
169     AliInfo("No Monte Carlo information in run");
170   }
171 }
172
173
174 //___________________________________________________________________________
175 void AliAnalysisTaskSingleMu::FinishTaskOutput()
176 {
177   //
178   /// Perform histogram normalization after the last analyzed event
179   /// but before merging
180   //
181
182   // Set the correct run limits for the histograms
183   // vs run number to cope with the merging
184   Int_t histoIndex = -1;
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   }
198 }
199
200
201 //___________________________________________________________________________
202 void AliAnalysisTaskSingleMu::UserCreateOutputObjects() 
203 {
204   //
205   /// Create output histograms
206   //
207   AliInfo(Form("   CreateOutputObjects of task %s\n", GetName()));
208
209   // initialize histogram lists
210   if ( ! fHistoList ) fHistoList = new TList();
211   if ( ! fHistoListMC ) fHistoListMC = new TList();
212   if ( ! fHistoListQA ) fHistoListQA = new TList();
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];
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
227   Int_t nPtBins = 60;
228   Double_t ptMin = 0., ptMax = 30.;
229   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
230
231   Int_t nRapidityBins = 25;
232   Double_t rapidityMin = -4.5, rapidityMax = -2.;
233   //TString rapidityName("Rapidity"), rapidityTitle("y"), rapidityUnits("");
234   TString rapidityName("Eta"), rapidityTitle("#eta"), rapidityUnits("");
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");
239   
240   Int_t nDcaBins = 30;
241   Double_t dcaMin = 0., dcaMax = 300.;
242   TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");
243
244   Int_t nVzBins = 40;
245   Double_t vzMin = -20., vzMax = 20.;
246   TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
247   
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("");
259   
260   Int_t nTrigClassBins = fTriggerClasses->GetEntries();
261   Double_t trigClassMin = 0.5, trigClassMax = (Double_t)nTrigClassBins + 0.5;
262   TString trigClassName("TrigClass"), trigClassTitle("Fired trigger class"), trigClassUnits("");
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("");
271
272   Int_t nPBins = 100;
273   Double_t pMin = 0., pMax = 400.;
274   TString pName("P"), pTitle("p"), pUnits("GeV/c");
275
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
295   // Multi-dimensional histo
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};
301
302   //TString stepTitle[kNsteps] = {"reconstructed", "in acceptance", "generated", "in acceptance (MC)"};
303   TString stepTitle[kNsteps] = {"reconstructed", "generated"};
304
305   // Create CF container
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);
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]);
320   }
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     }
335   }
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
345   /*
346   // MC kinematic cuts
347   AliCFTrackKineCuts *mcAccCuts = new AliCFTrackKineCuts();
348   mcAccCuts->SetNameTitle("mcAccCuts","MC-level acceptance cuts");
349   mcAccCuts->SetEtaRange(-4.,-2.5);
350   mcAccCuts->SetQAOn(fHistoListQA);
351
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);
357   */
358
359   TObjArray* mcGenList = new TObjArray(0) ;
360   mcGenList->AddLast(mcGenCuts);
361
362   /*
363   TObjArray* mcAccList = new TObjArray(0);
364   mcAccList->AddLast(mcAccCuts);
365
366   TObjArray* recAccList = new TObjArray(0);
367   recAccList->AddLast(recAccCuts);
368   */
369
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);    
384   }
385   //fCFManager->SetParticleCutsList(kStepAcceptance,recAccList);
386   fCFManager->SetParticleCutsList(kStepGeneratedMC,mcGenList);
387   //fCFManager->SetParticleCutsList(kStepAcceptanceMC,mcAccList);
388
389   // Summary histos
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
403   histoName = "histoMuonMultiplicity";
404   histoTitle = "Muon track multiplicity";
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";
414   histo2D = new TH2F(histoName.Data(), histoTitle.Data(), nTrigClassBins, trigClassMin, trigClassMax, nVzBins, vzMin, vzMax);
415   histo2D->GetXaxis()->SetTitle("Fired class");
416   SetAxisLabel(histo2D->GetXaxis());
417   histoTitle = Form("%s (%s)", vzTitle.Data(), vzUnits.Data());
418   histo2D->GetYaxis()->SetTitle(histoTitle.Data());
419   histoIndex = GetHistoIndex(kHistoEventVz);
420   fHistoList->AddAt(histo2D, histoIndex);
421
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
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()));
458       histoIndex = GetHistoIndex(kHistoPtResolutionMC, itrig, isrc);
459       fHistoListMC->AddAt(histo1D, histoIndex);
460     }
461   }
462
463   // Trees
464   if ( fFillTreeScaleDown > 0 ) {
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",
470                                        "Eta", "Rapidity", "Charge", "RAtAbsEnd",
471                                        "IPVx", "IPVy", "IPVz"};
472     TString leavesInt[kNvarInt] = {"MatchTrig", "IsMuon", "IsGhost", "LoCircuit", "PassPhysicsSelection", "NVtxContrib", "NspdTracklets", "IsPileupVertex"};
473     TString leavesChar[kNvarChar] = {"FiredTrigClass"};
474     TString leavesUInt[kNvarUInt] = {"BunchCrossNum", "OrbitNum", "PeriodNum", "RunNum"};
475     TString leavesFloatMC[kNvarFloatMC] = {"PxMC", "PyMC", "PzMC", "PtMC", "EtaMC", "RapidityMC", "VxMC", "VyMC", "VzMC"};
476     TString leavesIntMC[kNvarIntMC] = {"Pdg", "MotherType"};
477
478     TString treeName = GetOutputSlot(5)->GetContainer()->GetName();
479     Bool_t hasMC = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
480     TString treeTitle = ( hasMC ) ? "Single Mu" : "Single Mu MC";
481
482     OpenFile(5);
483     if ( ! fTreeSingleMu ) fTreeSingleMu = new TTree(treeName.Data(), treeTitle.Data());
484
485     for(Int_t itree=0; itree<2; itree++){
486       TParameter<Int_t>* par1 = new TParameter<Int_t>("fillTreeScaleDown",fFillTreeScaleDown);
487       TParameter<Int_t>* par2 = new TParameter<Int_t>("keepAllEvents",fKeepAll);
488       fTreeSingleMu->GetUserInfo()->Add(par1);
489       fTreeSingleMu->GetUserInfo()->Add(par2);
490       for(Int_t ivar=0; ivar<kNvarFloat; ivar++){
491         fTreeSingleMu->Branch(leavesFloat[ivar].Data(), &fVarFloat[ivar], Form("%s/F", leavesFloat[ivar].Data()));
492       }
493       for(Int_t ivar=0; ivar<kNvarInt; ivar++){
494         fTreeSingleMu->Branch(leavesInt[ivar].Data(), &fVarInt[ivar], Form("%s/I", leavesInt[ivar].Data()));
495       }
496       for(Int_t ivar=0; ivar<kNvarChar; ivar++){
497         TString addString = leavesChar[ivar] + "/C";
498         fTreeSingleMu->Branch(leavesChar[ivar].Data(), fVarChar[ivar], addString.Data());
499       }
500       for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
501         fTreeSingleMu->Branch(leavesUInt[ivar].Data(), &fVarUInt[ivar], Form("%s/i", leavesUInt[ivar].Data()));
502       }
503       if ( hasMC ){
504         for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
505           fTreeSingleMu->Branch(leavesFloatMC[ivar].Data(), &fVarFloatMC[ivar], Form("%s/F", leavesFloatMC[ivar].Data()));
506         }
507         for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
508           fTreeSingleMu->Branch(leavesIntMC[ivar].Data(), &fVarIntMC[ivar], Form("%s/I", leavesIntMC[ivar].Data()));
509         }
510       }
511     } // loop on trees
512     PostData(5, fTreeSingleMu);
513   } // fillNTuple
514   
515   PostData(1, fCFManager->GetParticleContainer());
516   PostData(2, fHistoList);
517   PostData(3, fHistoListQA);
518   PostData(4, fHistoListMC);
519
520 }
521
522 //________________________________________________________________________
523 void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/) 
524 {
525   //
526   /// Main loop
527   /// Called for each event
528   //
529
530   AliESDEvent* esdEvent = 0x0;
531   AliAODEvent* aodEvent = 0x0;
532
533   esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
534   if ( ! esdEvent ){
535     aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
536   }
537
538   if ( ! aodEvent && ! esdEvent ) {
539     AliError ("AOD or ESD event not found. Nothing done!");
540     return;
541   }
542
543   if ( ! fMCEvent && InputEvent()->GetEventType() != 7 ) return; // Run only on physics events!
544
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();
551   /*
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
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()));
559     }
560   }
561   */
562
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();
566   fVarInt[kVarIsPileup] = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8); // REMEMBER TO CHECK
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 ){
575     strncpy(fVarChar[kVarTrigMask], firedTrigClasses.Data(),255);
576
577     fVarInt[kVarPassPhysicsSelection] = ( esdEvent ) ? ((AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected() : ((AliAODInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
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;
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();
586
587     fVarFloat[kVarIPVx] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetX() : aodEvent->GetPrimaryVertex()->GetX();
588       fVarFloat[kVarIPVy] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetY() : aodEvent->GetPrimaryVertex()->GetY();
589   }
590
591   // Object declaration
592   AliMCParticle* mcPart = 0x0;
593   AliVParticle* track = 0x0;
594
595   Double_t containerInput[kNvars];
596   Int_t histoIndex = -1;
597
598   fCFManager->SetRecEventInfo(InputEvent());
599
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();
620       //containerInput[kHvarY]   = mcPart->Y();
621       containerInput[kHvarEta]   = mcPart->Eta();
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.;
628       containerInput[kHvarMotherType] = (Double_t)RecoTrackMother(mcPart);
629       containerInput[kHvarP] = mcPart->P();
630
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);
637         // if ( fCFManager->CheckParticleCuts(kStepAcceptanceMC,mcPart) ) fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptanceMC);
638       } // loop on trigger classes
639     } // loop on MC particles
640   } // is MC
641
642
643   Int_t trackLabel = -1;
644   Bool_t isGhost = kFALSE;
645   Int_t nGhosts = 0, nMuons = 0;
646
647   Int_t nTracks = ( esdEvent ) ? esdEvent->GetNumberOfMuonTracks() : aodEvent->GetNTracks();
648
649   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
650
651     // Get variables
652     if ( esdEvent ){
653       // ESD only info
654
655       track = esdEvent->GetMuonTrack(itrack);
656       isGhost = ( ((AliESDMuonTrack*)track)->ContainTriggerData() && ! ((AliESDMuonTrack*)track)->ContainTrackerData() );
657
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]);
667
668         fVarFloat[kVarXUncorrected] = ((AliESDMuonTrack*)track)->GetNonBendingCoorUncorrected();
669         fVarFloat[kVarYUncorrected] = ((AliESDMuonTrack*)track)->GetBendingCoorUncorrected();
670         fVarFloat[kVarZUncorrected] = ((AliESDMuonTrack*)track)->GetZUncorrected();
671
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();
682           fTreeSingleMu->Fill();
683         }
684         continue;
685       }
686     }
687     else {
688       track = aodEvent->GetTrack(itrack);
689       if ( ! ((AliAODTrack*)track)->IsMuonTrack() )
690         continue;
691     }
692
693     // Information for tracks in tracker
694     nMuons++;
695
696     fVarFloat[kVarPt] = track->Pt();
697     //fVarFloat[kVarRapidity] = track->Y();
698     fVarFloat[kVarEta] = track->Eta();
699     fVarFloat[kVarXatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetNonBendingCoorAtDCA() : ((AliAODTrack*)track)->XAtDCA();
700     fVarFloat[kVarYatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetBendingCoorAtDCA() : ((AliAODTrack*)track)->YAtDCA();
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();
706     fVarInt[kVarMatchTrig] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetMatchTrigger() : ((AliAODTrack*)track)->GetMatchTrigger();
707
708     if ( fillCurrentEventTree ){
709       fVarInt[kVarIsMuon] = nMuons;
710       fVarInt[kVarIsGhost] = 0;
711       //fVarFloat[kVarEta] = track->Eta();
712       fVarFloat[kVarRapidity] = track->Y();
713       fVarFloat[kVarPx] = track->Px();
714       fVarFloat[kVarPy] = track->Py();
715       fVarFloat[kVarPz] = track->Pz();
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();
719       fVarFloat[kVarPtAtDCA] = TMath::Sqrt(fVarFloat[kVarPxAtDCA]*fVarFloat[kVarPxAtDCA] + fVarFloat[kVarPyAtDCA]*fVarFloat[kVarPyAtDCA]);
720     }
721
722     fVarIntMC[kVarMotherType] = kUnknownPart;
723
724     // Monte Carlo part
725     if ( fMCEvent ) {
726
727       trackLabel = track->GetLabel();
728
729       AliMCParticle* matchedMCTrack = 0x0;
730
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         }
755       }
756       AliDebug(1, Form("Found mother %i", fVarIntMC[kVarMotherType]));
757     } // if use MC
758
759     containerInput[kHvarPt]  = fVarFloat[kVarPt];
760     //containerInput[kHvarY]   = fVarFloat[kVarRapidity];
761     containerInput[kHvarEta]   = fVarFloat[kVarEta];
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];
768     containerInput[kHvarIsGoodVtx] = (Double_t)isGoodVtxBin;
769     containerInput[kHvarMotherType] = (Double_t)fVarIntMC[kVarMotherType];
770     containerInput[kHvarP] = track->P();
771
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);
779       // if ( fCFManager->CheckParticleCuts(kStepAcceptance,track) ) fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptance);
780       ((TH2F*)fHistoList->At(histoIndex))->Fill(Form("%u",fVarUInt[kVarRunNumber]), containerInput[kHvarTrigClass], 1.);
781     } // loop on trigger classes
782
783     if ( fillCurrentEventTree ) fTreeSingleMu->Fill();
784   } // loop on tracks
785
786   if ( fillCurrentEventTree && fKeepAll &&  ( ( nMuons + nGhosts ) == 0 ) ) {
787     // Fill event also if there is not muon (when explicitely required)
788     fTreeSingleMu->Fill();
789   }
790
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     }
805
806     histoIndex = GetHistoIndex(kHistoMuonMultiplicity);
807     ((TH1F*)fHistoList->At(histoIndex))->Fill(nMuons, trigClassBin);
808     
809     if ( isGoodVtxBin == 1 ) {
810       histoIndex = GetHistoIndex(kHistoEventVz);
811       ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin,fVarFloat[kVarIPVz]);
812     }
813
814     histoIndex = GetHistoIndex(kHistoNeventsPerRun);
815     ((TH2F*)fHistoList->At(histoIndex))->Fill(Form("%u",fVarUInt[kVarRunNumber]), trigClassBin, 1.);
816   } // loop on trigger classes
817
818
819   // Post final data. It will be written to a file with option "RECREATE"
820   PostData(1, fCFManager->GetParticleContainer());
821   PostData(2, fHistoList);
822   PostData(3, fHistoListQA);
823   PostData(4, fHistoListMC);
824   if ( fFillTreeScaleDown > 0 )
825     PostData(5, fTreeSingleMu);
826 }
827
828 //________________________________________________________________________
829 void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
830   //
831   /// Draw some histogram at the end.
832   //
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() ) {
843     TString currName = GetName();
844     currName.Prepend("c1_");
845     TCanvas *c1_SingleMu = new TCanvas(currName.Data(),"Vz vs Pt",10,10,310,310);
846     c1_SingleMu->SetFillColor(10); c1_SingleMu->SetHighLightColor(10);
847     c1_SingleMu->SetLeftMargin(0.15); c1_SingleMu->SetBottomMargin(0.15);
848     TH2* histo = dynamic_cast<TH2*>(container->Project(kHvarPt,kHvarVz,kStepReconstructed));
849     currName = GetName();
850     currName.Prepend("hPtVz_");
851     histo->SetName(currName.Data());
852     histo->Draw("COLZ");
853   }
854
855 }
856
857 //________________________________________________________________________
858  void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
859                                                 Float_t var1, Float_t var2, Float_t var3)
860 {
861   //
862   /// Fill all histograms passing the trigger cuts
863   //
864
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   }
871
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   }
888 }
889
890 //________________________________________________________________________
891 Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcParticle)
892 {
893   //
894   /// Find track mother from kinematics
895   //
896
897   Int_t recoPdg = mcParticle->PdgCode();
898
899   // Track is not a muon
900   if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
901
902   Int_t imother = mcParticle->GetMother();
903
904   Bool_t isFirstMotherHF = kFALSE;
905   Int_t step = 0;
906
907   while ( imother >= 0 ) {
908     TParticle* part = ((AliMCParticle*)fMCEvent->GetTrack(imother))->Particle();
909
910     if ( part->GetUniqueID() == kPHadronic ) 
911       return kSecondaryMu;
912
913     Int_t absPdg = TMath::Abs(part->GetPdgCode());
914
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;
922
923     if ( isFirstMotherHF) {
924       if ( absPdg == 4 )
925         return kCharmMu;
926       else if ( absPdg == 5 )
927         return kBeautyMu;
928     }
929
930     imother = part->GetFirstMother();
931   }
932
933   return kPrimaryMu;
934 }
935
936 //________________________________________________________________________
937 Int_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;
945   }
946
947   return
948     kNsummaryHistosMC + 
949     kNtrackSources * kNtrigCuts * histoTypeIndex  + 
950     kNtrackSources * trigIndex  + 
951     srcIndex;
952 }
953
954 //________________________________________________________________________
955 Float_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
972
973 //________________________________________________________________________
974 void 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++){
985     fVarInt[ivar] = 0;
986   }
987   fVarInt[kVarMatchTrig] = -1;
988
989   if ( ! keepGlobal ){
990     for(Int_t ivar=0; ivar<kNvarChar; ivar++){
991       strncpy(fVarChar[ivar]," ",255);
992     }
993     for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
994       fVarUInt[ivar] = 0;
995     }
996   }
997   if ( fMCEvent ){
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 }
1006
1007 //________________________________________________________________________
1008 void AliAnalysisTaskSingleMu::SetTriggerClasses(TString triggerClasses)
1009 {
1010   /// Set trigger classes
1011   if ( fTriggerClasses )
1012     delete fTriggerClasses;
1013
1014   fTriggerClasses = triggerClasses.Tokenize(" ");
1015   fTriggerClasses->SetOwner();
1016
1017 }
1018
1019 //________________________________________________________________________
1020 void AliAnalysisTaskSingleMu::SetAxisLabel(TAxis* axis)
1021 {
1022   //
1023   /// Set the labels of trigger class axis
1024   //
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");
1030 }