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