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