PWG3/muon/AliAnalysisTaskSingleMu.cxx

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