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
  46 // STEER includes
  47 #include "AliLog.h"
  48
  49 #include "AliAODEvent.h"
  50 #include "AliAODTrack.h"
  51 #include "AliAODVertex.h"
  52
  53 #include "AliMCEvent.h"
  54 #include "AliMCParticle.h"
  55 #include "AliStack.h"
  56
  57 #include "AliESDInputHandler.h"
  58 #include "AliESDEvent.h"
  59 #include "AliESDMuonTrack.h"
  60 #include "AliAnalysisManager.h"
  61
  62 #include "AliAnalysisTaskSE.h"
  63 #include "AliAnalysisTaskSingleMu.h"
  64
  65 /// \cond CLASSIMP
  66 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
  67 /// \endcond
  68
  69
  70 //________________________________________________________________________
  71 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, Bool_t fillTree, Bool_t keepAll) :
  72   AliAnalysisTaskSE(name),
  73   fUseMC(0),
  74   fFillTree(fillTree),
  75   fKeepAll(keepAll),
  76   fHistoList(0),
  77   fHistoListMC(0),
  78   fTreeSingleMu(0),
  79   fTreeSingleMuMC(0),
  80   fVarFloat(0),
  81   fVarInt(0),
  82   fVarChar(0),
  83   fVarUInt(0),
  84   fVarFloatMC(0),
  85   fVarIntMC(0)
  86 {
  87   //
  88   /// Constructor.
  89   //
  90   if ( ! fFillTree )
  91     fKeepAll = kFALSE;
  92
  93   DefineOutput(1, TList::Class());
  94   DefineOutput(2, TList::Class());
  95
  96   if ( fFillTree )
  97     DefineOutput(3, TTree::Class());
  98 }
  99
 100
 101 //________________________________________________________________________
 102 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
 103 {
 104   //
 105   /// Destructor
 106   //
 107   delete fHistoList;
 108   delete fHistoListMC;
 109   delete fTreeSingleMu;
 110   delete fTreeSingleMuMC;
 111   delete fVarFloat;
 112   delete fVarInt;
 113   delete [] fVarChar;
 114   delete fVarUInt;
 115   delete fVarFloatMC;
 116   delete fVarIntMC;
 117 }
 118
 119
 120 //___________________________________________________________________________
 121 void AliAnalysisTaskSingleMu::NotifyRun()
 122 {
 123   //
 124   /// Use the event handler information to correctly fill the analysis flags:
 125   /// - check if Monte Carlo information is present
 126   //
 127
 128   if ( MCEvent() ) {
 129     fUseMC = kTRUE;
 130     AliInfo("Monte Carlo information is present");
 131   }
 132   else {
 133     AliInfo("No Monte Carlo information in run");
 134   }
 135 }
 136
 137
 138 //___________________________________________________________________________
 139 void AliAnalysisTaskSingleMu::UserCreateOutputObjects()
 140 {
 141   //
 142   /// Create output histograms
 143   //
 144   AliInfo(Form("   CreateOutputObjects of task %s\n", GetName()));
 145
 146   if ( fFillTree ) OpenFile(1);
 147
 148   // initialize histogram lists
 149   if(!fHistoList) fHistoList = new TList();
 150   if(!fHistoListMC) fHistoListMC = new TList();
 151
 152   // Init variables
 153   fVarFloat = new Float_t [kNvarFloat];
 154   fVarInt = new Int_t [kNvarInt];
 155   fVarChar = new Char_t *[kNvarChar];
 156   fVarUInt = new UInt_t [kNvarUInt];
 157   fVarFloatMC = new Float_t [kNvarFloatMC];
 158   fVarIntMC = new Int_t [kNvarIntMC];
 159
 160   Int_t nPtBins = 60;
 161   Float_t ptMin = 0., ptMax = 30.;
 162   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
 163
 164   Int_t nDcaBins = 100;
 165   Float_t dcaMin = 0., dcaMax = 200.;
 166   TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");
 167
 168   Int_t nVzBins = 60;
 169   Float_t vzMin = -30, vzMax = 30;
 170   TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
 171
 172   Int_t nEtaBins = 25;
 173   Float_t etaMin = -4.5, etaMax = -2.;
 174   TString etaName("Eta"), etaTitle("#eta"), etaUnits("a.u.");
 175
 176   Int_t nRapidityBins = 25;
 177   Float_t rapidityMin = -4.5, rapidityMax = -2.;
 178   TString rapidityName("Rapidity"), rapidityTitle("rapidity"), rapidityUnits("a.u.");
 179
 180   TString trigName[kNtrigCuts];
 181   trigName[kNoMatchTrig] = "NoMatch";
 182   trigName[kAllPtTrig]   = "AllPt";
 183   trigName[kLowPtTrig]   = "LowPt";
 184   trigName[kHighPtTrig]  = "HighPt";
 185
 186   TString srcName[kNtrackSources];
 187   srcName[kCharmMu]     = "Charm";
 188   srcName[kBeautyMu]    = "Beauty";
 189   srcName[kPrimaryMu]   = "Decay";
 190   srcName[kSecondaryMu] = "Secondary";
 191   srcName[kRecoHadron]  = "Hadrons";
 192   srcName[kUnknownPart] = "Unknown";
 193
 194   TH1F* histo1D = 0x0;
 195   TH2F* histo2D = 0x0;
 196   TString histoName, histoTitle;
 197   Int_t histoIndex = 0;
 198
 199   // 1D histos
 200   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 201     histoName = Form("%sDistrib%sTrig", ptName.Data(), trigName[itrig].Data());
 202     histoTitle = Form("%s distribution. Trigger: %s", ptTitle.Data(), trigName[itrig].Data());
 203     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nPtBins, ptMin, ptMax);
 204     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
 205     histoIndex = GetHistoIndex(kHistoPt, itrig);
 206     fHistoList->AddAt(histo1D, histoIndex);
 207   }
 208
 209   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 210     histoName = Form("%sDistrib%sTrig", dcaName.Data(), trigName[itrig].Data());
 211     histoTitle = Form("%s distribution. Trigger: %s", dcaTitle.Data(), trigName[itrig].Data());
 212     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nDcaBins, dcaMin, dcaMax);
 213     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
 214     histoIndex = GetHistoIndex(kHistoDCA, itrig);
 215     fHistoList->AddAt(histo1D, histoIndex);
 216   }
 217
 218   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 219     histoName = Form("%sDistrib%sTrig", vzName.Data(), trigName[itrig].Data());
 220     histoTitle = Form("%s distribution. Trigger: %s", vzTitle.Data(), trigName[itrig].Data());
 221     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nVzBins, vzMin, vzMax);
 222     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
 223     histoIndex = GetHistoIndex(kHistoVz, itrig);
 224     fHistoList->AddAt(histo1D, histoIndex);
 225   }
 226
 227   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 228     histoName = Form("%sDistrib%sTrig", etaName.Data(), trigName[itrig].Data());
 229     histoTitle = Form("%s distribution. Trigger: %s", etaTitle.Data(), trigName[itrig].Data());
 230     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nEtaBins, etaMin, etaMax);
 231     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", etaTitle.Data(), etaUnits.Data()));
 232     histoIndex = GetHistoIndex(kHistoEta, itrig);
 233     fHistoList->AddAt(histo1D, histoIndex);
 234   }
 235
 236   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 237     histoName = Form("%sDistrib%sTrig", rapidityName.Data(), trigName[itrig].Data());
 238     histoTitle = Form("%s distribution. Trigger: %s", rapidityTitle.Data(), trigName[itrig].Data());
 239     histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nRapidityBins, rapidityMin, rapidityMax);
 240     histo1D->GetXaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
 241     histoIndex = GetHistoIndex(kHistoRapidity, itrig);
 242     fHistoList->AddAt(histo1D, histoIndex);
 243   }
 244
 245   // 2D histos
 246   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 247     histoName = Form("%s%sDistrib%sTrig", dcaName.Data(), ptName.Data(), trigName[itrig].Data());
 248     histoTitle = Form("%s vs. %s distribution. Trigger: %s", dcaTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
 249     histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
 250                        nPtBins, ptMin, ptMax,
 251                        nDcaBins, dcaMin, dcaMax);
 252     histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
 253     histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
 254     histoIndex = GetHistoIndex(kHistoPtDCA, itrig);
 255     fHistoList->AddAt(histo2D, histoIndex);
 256   }
 257
 258   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 259     histoName = Form("%s%sDistrib%sTrig", vzName.Data(), ptName.Data(), trigName[itrig].Data());
 260     histoTitle = Form("%s vs. %s distribution. Trigger: %s", vzTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
 261     histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
 262                        nPtBins, ptMin, ptMax,
 263                        nVzBins, vzMin, vzMax);
 264     histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
 265     histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
 266     histoIndex = GetHistoIndex(kHistoPtVz, itrig);
 267     fHistoList->AddAt(histo2D, histoIndex);
 268   }
 269
 270   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 271     histoName = Form("%s%sDistrib%sTrig", rapidityName.Data(), ptName.Data(), trigName[itrig].Data());
 272     histoTitle = Form("%s vs. %s distribution. Trigger: %s", rapidityTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
 273     histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
 274                      nPtBins, ptMin, ptMax,
 275                      nRapidityBins, rapidityMin, rapidityMax);
 276     histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
 277     histo2D->GetYaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
 278     histoIndex = GetHistoIndex(kHistoPtRapidity, itrig);
 279     fHistoList->AddAt(histo2D, histoIndex);
 280   }
 281
 282   // Summary histos
 283   histoName = "histoMuonMultiplicity";
 284   histoTitle = "Muon track multiplicity";
 285   histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 10, -0.5, 10-0.5);
 286   //histo1D->GetXaxis()->SetBinLabel(1, "All events");
 287   //histo1D->GetXaxis()->SetBinLabel(2, "Muon events");
 288   histo1D->GetXaxis()->SetTitle("# of muons");
 289   histoIndex = GetHistoIndex(kNhistoTypes, 0) + kHistoMuonMultiplicity;
 290   fHistoList->AddAt(histo1D, histoIndex);
 291
 292   // MC histos
 293   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 294     for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
 295       histoName = Form("%sResolution%sTrig%s", ptName.Data(), trigName[itrig].Data(), srcName[isrc].Data());
 296       histoTitle = Form("%s resolution. Trigger: %s (%s)", ptTitle.Data(), trigName[itrig].Data(), srcName[isrc].Data());
 297       histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 100, -5., 5.);
 298       histo1D->GetXaxis()->SetTitle(Form("%s^{reco} - %s^{MC} (%s)", ptTitle.Data(), ptTitle.Data(), ptUnits.Data()));
 299       histoIndex = GetHistoIndex(kHistoPtResolutionMC, itrig, isrc);
 300       fHistoListMC->AddAt(histo1D, histoIndex);
 301     }
 302   }
 303
 304   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 305     for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
 306       histoName = Form("%s%sDistrib%sTrig%s", dcaName.Data(), ptName.Data(),
 307                        trigName[itrig].Data(), srcName[isrc].Data());
 308       histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", dcaTitle.Data(), ptTitle.Data(),
 309                         trigName[itrig].Data(), srcName[isrc].Data());
 310       histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
 311                          nPtBins, ptMin, ptMax,
 312                          nDcaBins, dcaMin, dcaMax);
 313       histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
 314       histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
 315       histoIndex = GetHistoIndex(kHistoPtDCAMC, itrig, isrc);
 316       fHistoListMC->AddAt(histo2D, histoIndex);
 317     }
 318   }
 319
 320   for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
 321     for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
 322       histoName = Form("%s%sDistrib%sTrig%s", vzName.Data(), ptName.Data(),
 323                        trigName[itrig].Data(), srcName[isrc].Data());
 324       histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", vzTitle.Data(), ptTitle.Data(),
 325                         trigName[itrig].Data(), srcName[isrc].Data());
 326       histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
 327                          nPtBins, ptMin, ptMax,
 328                          nVzBins, vzMin, vzMax);
 329       histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
 330       histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
 331       histoIndex = GetHistoIndex(kHistoPtVzMC, itrig, isrc);
 332       fHistoListMC->AddAt(histo2D, histoIndex);
 333     }
 334   }
 335
 336   // Trees
 337   if ( fFillTree ){
 338     TString leavesFloat[kNvarFloat] = {"Px", "Py", "Pz", "Pt",
 339                                        "PxAtDCA", "PyAtDCA", "PzAtDCA", "PtAtDCA",
 340                                        "PxUncorrected", "PyUncorrected", "PzUncorrected", "PtUncorrected",
 341                                        "XUncorrected", "YUncorrected", "ZUncorrected",
 342                                        "XatDCA", "YatDCA", "DCA",
 343                                        "Eta", "Rapidity", "Charge",
 344                                        "IPVx", "IPVy", "IPVz"};
 345     TString leavesInt[kNvarInt] = {"MatchTrig", "IsMuon", "IsGhost", "PassPhysicsSelection"};
 346     TString leavesChar[kNvarChar] = {"FiredTrigClass"};
 347     const Int_t charWidth[kNvarChar] = {255};
 348     TString leavesUInt[kNvarUInt] = {"BunchCrossNum", "OrbitNum", "PeriodNum", "RunNum"};
 349     TString leavesFloatMC[kNvarFloatMC] = {"PxMC", "PyMC", "PzMC", "PtMC", "EtaMC", "RapidityMC", "VxMC", "VyMC", "VzMC"};
 350     TString leavesIntMC[kNvarIntMC] = {"Pdg", "MotherType"};
 351
 352     if ( ! fTreeSingleMu ) fTreeSingleMu = new TTree("fTreeSingleMu", "Single Mu");
 353     if ( ! fTreeSingleMuMC ) fTreeSingleMuMC = new TTree("fTreeSingleMuMC", "Single Mu MC");
 354
 355     TTree* currTree[2] = {fTreeSingleMu, fTreeSingleMuMC};
 356     //TList* currList[2] = {fHistoList, fHistoListMC};
 357
 358     for(Int_t itree=0; itree<2; itree++){
 359       for(Int_t ivar=0; ivar<kNvarFloat; ivar++){
 360         currTree[itree]->Branch(leavesFloat[ivar].Data(), &fVarFloat[ivar], Form("%s/F", leavesFloat[ivar].Data()));
 361       }
 362       for(Int_t ivar=0; ivar<kNvarInt; ivar++){
 363         currTree[itree]->Branch(leavesInt[ivar].Data(), &fVarInt[ivar], Form("%s/I", leavesInt[ivar].Data()));
 364       }
 365       for(Int_t ivar=0; ivar<kNvarChar; ivar++){
 366         if ( itree == 0 ){
 367           fVarChar[ivar] = new Char_t [charWidth[ivar]];
 368         }
 369         TString addString = leavesChar[ivar] + "/C";
 370         currTree[itree]->Branch(leavesChar[ivar].Data(), fVarChar[ivar], addString.Data());
 371       }
 372       for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
 373         currTree[itree]->Branch(leavesUInt[ivar].Data(), &fVarUInt[ivar], Form("%s/i", leavesUInt[ivar].Data()));
 374       }
 375       if ( itree==1 ){
 376         for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
 377           currTree[itree]->Branch(leavesFloatMC[ivar].Data(), &fVarFloatMC[ivar], Form("%s/F", leavesFloatMC[ivar].Data()));
 378         }
 379         for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
 380           currTree[itree]->Branch(leavesIntMC[ivar].Data(), &fVarIntMC[ivar], Form("%s/I", leavesIntMC[ivar].Data()));
 381         }
 382       }
 383     } // loop on trees
 384   } // fillNTuple
 385 }
 386
 387 //________________________________________________________________________
 388 void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/)
 389 {
 390   //
 391   /// Main loop
 392   /// Called for each event
 393   //
 394
 395   AliESDEvent* esdEvent = 0x0;
 396   AliAODEvent* aodEvent = 0x0;
 397   AliMCEvent*  mcEvent  = 0x0;
 398
 399   esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
 400   if ( ! esdEvent ){
 401     fFillTree = kFALSE;
 402     aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
 403   }
 404   //else
 405   //aodEvent = AODEvent();
 406
 407   if ( ! aodEvent && ! esdEvent ) {
 408     AliError ("AOD or ESD event not found. Nothing done!");
 409     return;
 410   }
 411
 412   if ( ! fUseMC && InputEvent()->GetEventType() !=7 ) return; // Run only on physics events!
 413
 414   if ( fFillTree ){
 415     Reset(kFALSE);
 416     fVarInt[kVarPassPhysicsSelection] = ((AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
 417     strcpy(fVarChar[kVarTrigMask], esdEvent->GetFiredTriggerClasses().Data());
 418
 419     // Small workaround: in MC the bunch ID are not properly set and the timestamp is in seconds
 420     // So fill bunchCrossing with the read timestamp
 421     //    fill the orbit and period number with a timestamp created while reading the run
 422     TTimeStamp ts;
 423     fVarUInt[kVarBunchCrossNumber] = ( fUseMC ) ? (UInt_t)Entry() : esdEvent->GetBunchCrossNumber();
 424     fVarUInt[kVarOrbitNumber] = ( fUseMC ) ? ts.GetNanoSec() : esdEvent->GetOrbitNumber();
 425     fVarUInt[kVarPeriodNumber] = ( fUseMC ) ? ts.GetTime() : esdEvent->GetPeriodNumber();
 426     fVarUInt[kVarRunNumber] = esdEvent->GetRunNumber();
 427
 428     fVarFloat[kVarIPVx] = esdEvent->GetPrimaryVertex()->GetX();
 429     fVarFloat[kVarIPVy] = esdEvent->GetPrimaryVertex()->GetY();
 430     fVarFloat[kVarIPVz] = esdEvent->GetPrimaryVertex()->GetZ();
 431   }
 432
 433   if ( fUseMC ) mcEvent = MCEvent();
 434
 435   // Object declaration
 436   AliMCParticle* mcTrack = 0x0;
 437
 438   Int_t trackLabel = -1;
 439
 440   Int_t nTracks = ( esdEvent ) ? esdEvent->GetNumberOfMuonTracks() : aodEvent->GetNTracks();
 441
 442   Bool_t isGhost = kFALSE;
 443   Int_t nGhosts = 0, nMuons = 0;
 444
 445   AliVParticle* track = 0x0;
 446
 447   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
 448
 449     // Get variables
 450     if ( esdEvent ){
 451       if (itrack>0) Reset(kTRUE);
 452
 453       track = esdEvent->GetMuonTrack(itrack);
 454       isGhost = ( ((AliESDMuonTrack*)track)->ContainTriggerData() && ! ((AliESDMuonTrack*)track)->ContainTrackerData() );
 455
 456       // ESD only info
 457       fVarFloat[kVarPxUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaXUncorrected()) : ((AliESDMuonTrack*)track)->PxUncorrected();
 458       fVarFloat[kVarPyUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaYUncorrected()) : ((AliESDMuonTrack*)track)->PyUncorrected();
 459       fVarFloat[kVarPzUncorrected] = ( isGhost ) ? -1 : ((AliESDMuonTrack*)track)->PzUncorrected();
 460       fVarFloat[kVarPtUncorrected] = TMath::Sqrt(
 461         fVarFloat[kVarPxUncorrected] * fVarFloat[kVarPxUncorrected] +
 462         fVarFloat[kVarPyUncorrected] * fVarFloat[kVarPyUncorrected]);
 463
 464       fVarFloat[kVarXUncorrected] = ((AliESDMuonTrack*)track)->GetNonBendingCoorUncorrected();
 465       fVarFloat[kVarYUncorrected] = ((AliESDMuonTrack*)track)->GetBendingCoorUncorrected();
 466       fVarFloat[kVarZUncorrected] = ((AliESDMuonTrack*)track)->GetZUncorrected();
 467
 468       fVarInt[kVarMatchTrig] = ((AliESDMuonTrack*)track)->GetMatchTrigger();
 469
 470       // If is ghost fill only a partial information
 471       if ( isGhost ){
 472         fVarInt[kVarIsMuon] = 0;
 473         nGhosts++;
 474         fVarInt[kVarIsGhost] = nGhosts;
 475
 476         if ( ! fUseMC ) fTreeSingleMu->Fill();
 477         else fTreeSingleMuMC->Fill();
 478
 479         continue;
 480       }
 481     }
 482     else {
 483       track = aodEvent->GetTrack(itrack);
 484       if ( ! ((AliAODTrack*)track)->IsMuonTrack() )
 485         continue;
 486
 487       //Reset(kTRUE);
 488     }
 489
 490     // Information for tracks in tracker
 491     nMuons++;
 492     fVarInt[kVarIsMuon] = nMuons;
 493     fVarInt[kVarIsGhost] = 0;
 494
 495     fVarFloat[kVarPt] = track->Pt();
 496     fVarFloat[kVarXatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetNonBendingCoorAtDCA() : ((AliAODTrack*)track)->XAtDCA();
 497     fVarFloat[kVarYatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetBendingCoorAtDCA() : ((AliAODTrack*)track)->YAtDCA();
 498     fVarFloat[kVarDCA] = TMath::Sqrt( fVarFloat[kVarXatDCA] * fVarFloat[kVarXatDCA] + fVarFloat[kVarYatDCA] * fVarFloat[kVarYatDCA] );
 499     fVarFloat[kVarEta] = track->Eta();
 500     fVarFloat[kVarRapidity] = track->Y();
 501     trackLabel = track->GetLabel();
 502     fVarInt[kVarMatchTrig] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetMatchTrigger() : ((AliAODTrack*)track)->GetMatchTrigger();
 503
 504     // Fill histograms
 505     FillTriggerHistos(kHistoPt,       fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt]);
 506     FillTriggerHistos(kHistoDCA,      fVarInt[kVarMatchTrig], -1, fVarFloat[kVarDCA]);
 507     FillTriggerHistos(kHistoVz,       fVarInt[kVarMatchTrig], -1, fVarFloat[kVarIPVz]);
 508     FillTriggerHistos(kHistoEta,      fVarInt[kVarMatchTrig], -1, fVarFloat[kVarEta]);
 509     FillTriggerHistos(kHistoRapidity, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarRapidity]);
 510
 511     FillTriggerHistos(kHistoPtDCA,      fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt], fVarFloat[kVarDCA]);
 512     FillTriggerHistos(kHistoPtVz,       fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt], fVarFloat[kVarIPVz]);
 513     FillTriggerHistos(kHistoPtRapidity, fVarInt[kVarMatchTrig], -1, fVarFloat[kVarPt], fVarFloat[kVarRapidity]);
 514
 515     if ( fFillTree ){
 516       fVarFloat[kVarPx] = track->Px();
 517       fVarFloat[kVarPy] = track->Py();
 518       fVarFloat[kVarPz] = track->Pz();
 519       fVarFloat[kVarPxAtDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->PxAtDCA() : ((AliAODTrack*)track)->PxAtDCA();
 520       fVarFloat[kVarPyAtDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->PyAtDCA() : ((AliAODTrack*)track)->PyAtDCA();
 521       fVarFloat[kVarPzAtDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->PzAtDCA() : ((AliAODTrack*)track)->PzAtDCA();
 522       fVarFloat[kVarPtAtDCA] = TMath::Sqrt(fVarFloat[kVarPxAtDCA]*fVarFloat[kVarPxAtDCA] + fVarFloat[kVarPyAtDCA]*fVarFloat[kVarPyAtDCA]);
 523
 524       fVarFloat[kVarCharge] = track->Charge();
 525
 526       if ( ! fUseMC ) fTreeSingleMu->Fill();
 527     }
 528
 529     // Monte Carlo part
 530     if ( ! fUseMC ) continue;
 531
 532     fVarIntMC[kVarMotherType] = kUnknownPart;
 533
 534     AliMCParticle* matchedMCTrack = 0x0;
 535
 536     Int_t nMCtracks = mcEvent->GetNumberOfTracks();
 537     for(Int_t imctrack=0; imctrack<nMCtracks; imctrack++){
 538       mcTrack = (AliMCParticle*)mcEvent->GetTrack(imctrack);
 539       if ( trackLabel == mcTrack->GetLabel() ) {
 540         matchedMCTrack = mcTrack;
 541         break;
 542       }
 543     } // loop on MC tracks
 544
 545     if ( matchedMCTrack )
 546       fVarIntMC[kVarMotherType] = RecoTrackMother(matchedMCTrack, mcEvent);
 547
 548     AliDebug(1, Form("Found mother %i", fVarIntMC[kVarMotherType]));
 549
 550
 551     FillTriggerHistos(kHistoPtDCAMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt], fVarFloat[kVarDCA]);
 552     FillTriggerHistos(kHistoPtVzMC,  fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt], fVarFloat[kVarIPVz]);
 553
 554     if ( matchedMCTrack ) {
 555       fVarFloatMC[kVarPtMC] = matchedMCTrack->Pt();
 556       FillTriggerHistos(kHistoPtResolutionMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt] - fVarFloatMC[kVarPtMC]);
 557       if ( fFillTree ){
 558         fVarFloatMC[kVarPxMC] = matchedMCTrack->Px();
 559         fVarFloatMC[kVarPyMC] = matchedMCTrack->Py();
 560         fVarFloatMC[kVarPzMC] = matchedMCTrack->Pz();
 561         fVarFloatMC[kVarEtaMC] = matchedMCTrack->Eta();
 562         fVarFloatMC[kVarRapidityMC] = matchedMCTrack->Y();
 563         fVarFloatMC[kVarVxMC] = matchedMCTrack->Xv();
 564         fVarFloatMC[kVarVyMC] = matchedMCTrack->Yv();
 565         fVarFloatMC[kVarVzMC] = matchedMCTrack->Zv();
 566         fVarIntMC[kVarPdg] = matchedMCTrack->PdgCode();
 567       }
 568     }
 569     if ( fFillTree ) fTreeSingleMuMC->Fill();
 570   } // loop on tracks
 571
 572   if ( fKeepAll &&  ( ( fVarInt[kVarIsMuon] + fVarInt[kVarIsGhost] ) == 0 ) ) {
 573     // Fill event also if there is not muon (when explicitely required)
 574     if ( ! fUseMC ) fTreeSingleMu->Fill();
 575     else fTreeSingleMuMC->Fill();
 576   }
 577
 578   if( strstr(fVarChar[kVarTrigMask],"MUON") && fVarInt[kVarIsMuon]==0 )
 579     printf("WARNING: Muon trigger does not match tracker!\n");
 580
 581   Int_t histoIndex = GetHistoIndex(kNhistoTypes, 0) + kHistoMuonMultiplicity;
 582   ((TH1F*)fHistoList->At(histoIndex))->Fill(fVarInt[kVarIsMuon]);
 583   //if ( fVarInt[kVarIsMuon]>0 ) ((TH1F*)fHistoList->At(histoIndex))->Fill(2);
 584
 585   // Post final data. It will be written to a file with option "RECREATE"
 586   PostData(1, fHistoList);
 587   if ( fUseMC ) PostData(2, fHistoListMC);
 588   if ( fFillTree ){
 589     if ( fUseMC )
 590       PostData(3, fTreeSingleMuMC);
 591     else
 592       PostData(3, fTreeSingleMu);
 593   }
 594 }
 595
 596 //________________________________________________________________________
 597 void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
 598   //
 599   /// Draw some histogram at the end.
 600   //
 601   if (!gROOT->IsBatch()) {
 602     fHistoList = dynamic_cast<TList*> (GetOutputData(1));
 603     TCanvas *c1_SingleMu = new TCanvas("c1_SingleMu","Vz vs Pt",10,10,310,310);
 604     c1_SingleMu->SetFillColor(10); c1_SingleMu->SetHighLightColor(10);
 605     c1_SingleMu->SetLeftMargin(0.15); c1_SingleMu->SetBottomMargin(0.15);
 606     Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
 607     ((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
 608   }
 609 }
 610
 611 //________________________________________________________________________
 612  void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
 613                                                 Float_t var1, Float_t var2, Float_t var3)
 614 {
 615   //
 616   /// Fill all histograms passing the trigger cuts
 617   //
 618
 619   Int_t nTrigs = TMath::Max(1, matchTrig);
 620   TArrayI trigToFill(nTrigs);
 621   trigToFill[0] = matchTrig;
 622   for(Int_t itrig = 1; itrig < matchTrig; itrig++){
 623     trigToFill[itrig] = itrig;
 624   }
 625
 626   TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;
 627
 628   TString className;
 629   for(Int_t itrig = 0; itrig < nTrigs; itrig++){
 630     Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
 631     className = histoList->At(currIndex)->ClassName();
 632     AliDebug(3, Form("matchTrig %i  Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
 633     if (className.Contains("1"))
 634       ((TH1F*)histoList->At(currIndex))->Fill(var1);
 635     else if (className.Contains("2"))
 636       ((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
 637     else if (className.Contains("3"))
 638       ((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
 639     else
 640       AliWarning("Histogram not filled");
 641   }
 642 }
 643
 644 //________________________________________________________________________
 645 Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
 646 {
 647   //
 648   /// Find track mother from kinematics
 649   //
 650
 651   Int_t recoPdg = mcTrack->PdgCode();
 652
 653   Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;
 654
 655   // Track is not a muon
 656   if ( ! isMuon ) return kRecoHadron;
 657
 658   Int_t imother = mcTrack->GetMother();
 659
 660   if ( imother<0 )
 661     return kPrimaryMu; // Drell-Yan Muon
 662
 663   Int_t igrandma = imother;
 664
 665   AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
 666   Int_t motherPdg = motherPart->PdgCode();
 667
 668   // Track is an heavy flavor muon
 669   Int_t absPdg = TMath::Abs(motherPdg);
 670   if(absPdg/100==5 || absPdg/1000==5) {
 671     return kBeautyMu;
 672   }
 673   if(absPdg/100==4 || absPdg/1000==4){
 674     Int_t newMother = -1;
 675     igrandma = imother;
 676     AliInfo("\nFound candidate B->C->mu. History:\n");
 677     mcTrack->Print();
 678     printf("- %2i   ", igrandma);
 679     motherPart->Print();
 680     Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
 681     while ( absGrandMotherPdg > 10 ) {
 682       igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
 683       if( igrandma < 0 ) break;
 684       printf("- %2i   ", igrandma);
 685       mcEvent->GetTrack(igrandma)->Print();
 686       absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
 687     }
 688
 689     if (absGrandMotherPdg==5) newMother = kBeautyMu; // Charm from beauty
 690     else if (absGrandMotherPdg==4) newMother = kCharmMu;
 691
 692     if(newMother<0) {
 693       AliWarning("Mother not correctly found! Set to charm!\n");
 694       newMother = kCharmMu;
 695     }
 696
 697     return newMother;
 698   }
 699
 700   Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
 701
 702   // Track is a bkg. muon
 703   if (imother<nPrimaries) {
 704     return kPrimaryMu;
 705   }
 706   else {
 707     return kSecondaryMu;
 708   }
 709 }
 710
 711 //________________________________________________________________________
 712 Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
 713  {
 714    //
 715    /// Get histogram index in the list
 716    //
 717    if ( srcIndex < 0 ) return kNtrigCuts * histoTypeIndex + trigIndex;
 718
 719    return
 720      kNtrackSources * kNtrigCuts * histoTypeIndex  +
 721      kNtrackSources * trigIndex  +
 722      srcIndex;
 723 }
 724
 725 //________________________________________________________________________
 726 void AliAnalysisTaskSingleMu::Reset(Bool_t keepGlobal)
 727 {
 728   //
 729   /// Reset variables
 730   //
 731   Int_t lastVarFloat = ( keepGlobal ) ? kVarIPVx : kNvarFloat;
 732   for(Int_t ivar=0; ivar<lastVarFloat; ivar++){
 733     fVarFloat[ivar] = 0.;
 734   }
 735   Int_t lastVarInt = ( keepGlobal ) ? kVarPassPhysicsSelection : kNvarInt;
 736   for(Int_t ivar=0; ivar<lastVarInt; ivar++){
 737     fVarInt[ivar] = -1;
 738   }
 739   fVarInt[kVarIsMuon] = 0;
 740   fVarInt[kVarIsGhost] = 0;
 741
 742   if ( ! keepGlobal ){
 743     for(Int_t ivar=0; ivar<kNvarChar; ivar++){
 744       //sprintf(fVarChar[ivar],"%253s","");
 745       sprintf(fVarChar[ivar]," ");
 746     }
 747     for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
 748       fVarUInt[ivar] = 0;
 749     }
 750   }
 751   if ( fUseMC ){
 752     for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
 753       fVarFloatMC[ivar] = 0.;
 754     }
 755     for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
 756       fVarIntMC[ivar] = -1;
 757     }
 758   }
 759 }