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