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 }