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 and CF containers.
  22 /// The macro class can run on AODs or ESDs.
  23 /// If Monte Carlo information is present, some basics checks are performed.
  24 ///
  25 /// \author Diego Stocco
  26 //-----------------------------------------------------------------------------
  27
  28 #define AliAnalysisTaskSingleMu_cxx
  29
  30 #include "AliAnalysisTaskSingleMu.h"
  31
  32 // ROOT includes
  33 #include "TROOT.h"
  34 #include "TH1.h"
  35 #include "TH2.h"
  36 #include "TAxis.h"
  37 #include "TCanvas.h"
  38 #include "TLegend.h"
  39 #include "TMath.h"
  40 #include "TObjString.h"
  41 #include "TObjArray.h"
  42 #include "TF1.h"
  43 #include "TStyle.h"
  44 //#include "TMCProcess.h"
  45 #include "TArrayI.h"
  46 #include "TPaveStats.h"
  47 #include "TFitResultPtr.h"
  48
  49 // STEER includes
  50 #include "AliAODEvent.h"
  51 #include "AliAODTrack.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMCEvent.h"
  54 #include "AliMCParticle.h"
  55 #include "AliESDEvent.h"
  56 #include "AliESDMuonTrack.h"
  57 #include "AliVHeader.h"
  58 #include "AliAODMCHeader.h"
  59 #include "AliStack.h"
  60
  61 // ANALYSIS includes
  62 #include "AliAnalysisManager.h"
  63
  64 // CORRFW includes
  65 #include "AliCFContainer.h"
  66 #include "AliCFGridSparse.h"
  67 #include "AliCFEffGrid.h"
  68
  69 // PWG includes
  70 #include "AliVAnalysisMuon.h"
  71 #include "AliMergeableCollection.h"
  72 #include "AliCounterCollection.h"
  73 #include "AliMuonEventCuts.h"
  74 #include "AliMuonTrackCuts.h"
  75 #include "AliAnalysisMuonUtility.h"
  76
  77
  78 /// \cond CLASSIMP
  79 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
  80 /// \endcond
  81
  82
  83 //________________________________________________________________________
  84 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu() :
  85   AliVAnalysisMuon(),
  86   fThetaAbsKeys(0x0)
  87 {
  88   /// Default ctor.
  89 }
  90
  91 //________________________________________________________________________
  92 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
  93   AliVAnalysisMuon(name, cuts),
  94   fThetaAbsKeys(0x0)
  95 {
  96   //
  97   /// Constructor.
  98   //
  99   TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
 100   fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
 101 }
 102
 103
 104 //________________________________________________________________________
 105 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
 106 {
 107   //
 108   /// Destructor
 109   //
 110
 111   delete fThetaAbsKeys;
 112 }
 113
 114 //___________________________________________________________________________
 115 void AliAnalysisTaskSingleMu::MyUserCreateOutputObjects()
 116 {
 117
 118   TH1* histo = 0x0;
 119   TString histoName = "", histoTitle = "";
 120
 121   Int_t nVzBins = 40;
 122   Double_t vzMin = -20., vzMax = 20.;
 123   TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
 124
 125   histoName = "hIpVtx";
 126   histo = new TH1F(histoName.Data(), histoName.Data(), nVzBins, vzMin, vzMax);
 127   histo->SetXTitle("v_{z} (cm)");
 128   AddObjectToCollection(histo, kIPVz);
 129
 130
 131   Int_t nPtBins = 80;
 132   Double_t ptMin = 0., ptMax = 80.;
 133   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
 134
 135   Int_t nEtaBins = 25;
 136   Double_t etaMin = -4.5, etaMax = -2.;
 137   TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
 138
 139   Int_t nPhiBins = 36;
 140   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
 141   TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
 142
 143   Int_t nChargeBins = 2;
 144   Double_t chargeMin = -2, chargeMax = 2.;
 145   TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
 146
 147   Int_t nThetaAbsEndBins = 2;
 148   Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 1.5;
 149   TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");
 150
 151   Int_t nMotherTypeBins = kNtrackSources;
 152   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
 153   TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
 154
 155   Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nVzBins, nChargeBins, nThetaAbsEndBins, nMotherTypeBins};
 156   Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, vzMin, chargeMin, thetaAbsEndMin, motherTypeMin};
 157   Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, vzMax, chargeMax, thetaAbsEndMax, motherTypeMax};
 158   TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, vzTitle, chargeTitle, thetaAbsEndTitle, motherTypeTitle};
 159   TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, vzUnits, chargeUnits, thetaAbsEndUnits, motherTypeUnits};
 160
 161   AliCFContainer* cfContainer = new AliCFContainer("SingleMuContainer","Container for tracks",kNsteps,kNvars,nbins);
 162
 163   for ( Int_t idim = 0; idim<kNvars; idim++){
 164     histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
 165     histoTitle.ReplaceAll("()","");
 166
 167     cfContainer->SetVarTitle(idim, histoTitle.Data());
 168     cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
 169   }
 170
 171   TString stepTitle[kNsteps] = {"reconstructed", "generated"};
 172
 173   TAxis* currAxis = 0x0;
 174   for (Int_t istep=0; istep<kNsteps; istep++){
 175     cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
 176     AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
 177
 178     currAxis = gridSparse->GetAxis(kHvarMotherType);
 179     for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
 180       currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
 181     }
 182   }
 183
 184   AddObjectToCollection(cfContainer, kTrackContainer);
 185
 186   fMuonTrackCuts->Print("mask");
 187 }
 188
 189 //________________________________________________________________________
 190 void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
 191 {
 192   //
 193   /// Fill output objects
 194   //
 195
 196   Double_t ipVz = AliAnalysisMuonUtility::GetVertexSPD(InputEvent())->GetZ();
 197   Double_t ipVzMC = AliAnalysisMuonUtility::IsMCEvent(InputEvent(),MCEvent()) ? AliAnalysisMuonUtility::GetMCVertexZ(InputEvent(),MCEvent()) : 0.;
 198
 199   for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
 200     TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
 201     ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz);
 202   }
 203
 204   Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes
 205   if ( isPileupFromSPD ) return;
 206
 207   Double_t containerInput[kNvars];
 208   AliVParticle* track = 0x0;
 209
 210   for ( Int_t istep = 0; istep<2; ++istep ) {
 211     Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : AliAnalysisMuonUtility::GetNMCTracks(InputEvent(),MCEvent());
 212     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
 213       track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : AliAnalysisMuonUtility::GetMCTrack(itrack,InputEvent(),MCEvent());
 214
 215       Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
 216       if ( ! isSelected ) continue;
 217
 218       // In W simulations with Pythia, sometimes muon is stored twice.
 219       // Remove muon in case it has another muon as daugther
 220       if ( istep == kStepGeneratedMC ) {
 221         Int_t firstDaughter = AliAnalysisMuonUtility::GetDaughterIndex(track, 0);
 222         if ( firstDaughter >= 0 ) {
 223           Bool_t hasMuonDaughter = kFALSE;
 224           Int_t lastDaughter = AliAnalysisMuonUtility::GetDaughterIndex(track, 1);
 225           for ( Int_t idaugh=firstDaughter; idaugh<=lastDaughter; idaugh++ ) {
 226             AliVParticle* currTrack = AliAnalysisMuonUtility::GetMCTrack(idaugh,InputEvent(),MCEvent());
 227             if ( currTrack->PdgCode() == track->PdgCode() ) {
 228               hasMuonDaughter = kTRUE;
 229               break;
 230             }
 231           }
 232           if ( hasMuonDaughter ) {
 233             AliDebug(1, Form("Current muon (%i) has muon daughter: rejecting it", itrack));
 234             continue;
 235           }
 236         }
 237       }
 238
 239       Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
 240
 241       Double_t thetaAbsEndDeg = 0;
 242       if ( istep == kStepReconstructed ) {
 243         thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
 244       }
 245       else {
 246         thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg();
 247       }
 248       Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
 249
 250       containerInput[kHvarPt]         = track->Pt();
 251       containerInput[kHvarEta]        = track->Eta();
 252       containerInput[kHvarPhi]        = track->Phi();
 253       containerInput[kHvarVz]         = ( istep == kStepReconstructed ) ? ipVz : ipVzMC;
 254       containerInput[kHvarCharge]     = track->Charge()/3.;
 255       containerInput[kHvarThetaAbs]   = (Double_t)thetaAbsBin;
 256       containerInput[kHvarMotherType] = (Double_t)trackSrc;
 257
 258       for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
 259         TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
 260         if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
 261         ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep);
 262       } // loop on selected trigger classes
 263     } // loop on tracks
 264   } // loop on container steps
 265 }
 266
 267
 268 //________________________________________________________________________
 269 void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
 270   //
 271   /// Draw some histograms at the end.
 272   //
 273
 274   AliVAnalysisMuon::Terminate("");
 275
 276   if ( ! fMergeableCollection ) return;
 277
 278   TString physSel = fTerminateOptions->At(0)->GetName();
 279   TString trigClassName = fTerminateOptions->At(1)->GetName();
 280   TString centralityRange = fTerminateOptions->At(2)->GetName();
 281   TString furtherOpt = fTerminateOptions->At(3)->GetName();
 282
 283   TString minBiasTrig = "";
 284   TObjArray* optArr = furtherOpt.Tokenize(" ");
 285   TString currName = "";
 286   for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
 287     currName = optArr->At(iopt)->GetName();
 288     if ( currName.Contains("-B-") || currName.Contains("ANY") ) minBiasTrig = currName;
 289   }
 290   delete optArr;
 291
 292   furtherOpt.ToUpper();
 293   Bool_t plotChargeAsymmetry = furtherOpt.Contains("ASYM");
 294
 295   AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,centralityRange,"SingleMuContainer") );
 296   if ( ! cfContainer ) return;
 297
 298   AliCFEffGrid* effSparse = new AliCFEffGrid(Form("eff%s", cfContainer->GetName()),Form("Efficiency %s", cfContainer->GetTitle()),*cfContainer);
 299   effSparse->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC);
 300
 301   AliCFGridSparse* gridSparseArray[3] = {effSparse->GetNum(), effSparse->GetDen(), effSparse};
 302   TString gridSparseName[3] = {cfContainer->GetStepTitle(kStepReconstructed), cfContainer->GetStepTitle(kStepGeneratedMC), "Efficiency"};
 303
 304   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
 305 //  TString allSrcNames = "";
 306 //  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 307 //    if ( ! allSrcNames.IsNull() ) allSrcNames.Append(" ");
 308 //    allSrcNames += fSrcKeys->At(isrc)->GetName();
 309 //  }
 310
 311   TCanvas* can = 0x0;
 312   Int_t xshift = 100;
 313   Int_t yshift = 100;
 314   Int_t igroup1 = -1;
 315   Int_t igroup2 = 0;
 316
 317   Bool_t isMC = furtherOpt.Contains("MC");
 318   Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
 319   Int_t lastSrc  = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
 320   if ( ! isMC ) srcColors[kUnidentified] = 1;
 321
 322   TString histoName = "", histoPattern = "", drawOpt = "";
 323   ////////////////
 324   // Kinematics //
 325   ////////////////
 326   TCanvas* canKine[3] = {0x0, 0x0, 0x0};
 327   TLegend* legKine[3] = {0x0, 0x0, 0x0};
 328   for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
 329     for ( Int_t icharge=0; icharge<2; ++icharge ) {
 330       for ( Int_t igrid=0; igrid<3; ++igrid ) {
 331         if ( gridSparseArray[igrid]->GetEntries() == 0. ) break;
 332         if ( gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) {
 333           SetSparseRange(gridSparseArray[igrid], kHvarEta, "", -3.999, -2.501);
 334           SetSparseRange(gridSparseArray[igrid], kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
 335           SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", icharge+1, icharge+1, "USEBIN");
 336         }
 337         if ( ! canKine[igrid] ) {
 338           igroup1++;
 339           igroup2 = 0;
 340           currName = Form("%s_proj_%s", GetName(), gridSparseName[igrid].Data());
 341           canKine[igrid] = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 342           canKine[igrid]->Divide(2,2);
 343           legKine[igrid] = new TLegend(0.6, 0.6, 0.8, 0.8);
 344           igroup2++;
 345         }
 346         for ( Int_t iproj=0; iproj<4; ++iproj ) {
 347           canKine[igrid]->cd(iproj+1);
 348           if ( ( iproj == kHvarPt || iproj == kHvarVz ) && gridSparseArray[igrid]->IsA() != AliCFEffGrid::Class() ) gPad->SetLogy();
 349           TH1* projHisto = gridSparseArray[igrid]->Project(iproj);
 350           projHisto->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), fSrcKeys->At(isrc)->GetName(), fChargeKeys->At(icharge)->GetName()));
 351           if ( projHisto->GetEntries() == 0 ) continue;
 352           Bool_t isFirst = ( gPad->GetListOfPrimitives()->GetEntries() == 0 );
 353           drawOpt = isFirst ? "e" : "esames";
 354           //if ( isrc == kUnidentified && ! drawOpt.Contains("same") ) isMC = kFALSE;
 355           //if ( ! isMC ) srcColors[kUnidentified] = 1;
 356           projHisto->SetLineColor(srcColors[isrc]);
 357           projHisto->SetMarkerColor(srcColors[isrc]);
 358           projHisto->SetMarkerStyle(20+4*icharge);
 359           projHisto->Draw(drawOpt.Data());
 360           gPad->Update();
 361           TPaveStats* paveStats = (TPaveStats*)projHisto->FindObject("stats");
 362           if ( paveStats ) paveStats->SetTextColor(srcColors[isrc]);
 363           if ( iproj == 0 ) {
 364             TString legEntry = fChargeKeys->At(icharge)->GetName();
 365             if ( isMC ) legEntry += Form(" %s", fSrcKeys->At(isrc)->GetName());
 366             legKine[igrid]->AddEntry(projHisto,legEntry.Data(), "lp");
 367           }
 368         } // loop on grid sparse
 369       } // loop on projections
 370     } // loop on mu charge
 371   } // loop on track sources
 372
 373
 374   for ( Int_t igrid=0; igrid<3; igrid++ ) {
 375     if ( ! canKine[igrid] ) continue;
 376     canKine[igrid]->cd(1);
 377     legKine[igrid]->Draw("same");
 378     if ( gridSparseArray[igrid]->IsA() == AliCFEffGrid::Class() ) continue;
 379     SetSparseRange(gridSparseArray[igrid], kHvarCharge, "", 1, gridSparseArray[igrid]->GetAxis(kHvarCharge)->GetNbins(), "USEBIN"); // Reset range
 380   } // loop on container steps
 381
 382   // Plot charge asymmetry or mu+/mu-
 383   TString basePlotName = plotChargeAsymmetry ? "ChargeAsymmetry" : "ChargeRatio";
 384   for ( Int_t igrid=0; igrid<2; igrid++ ) {
 385     if ( ! canKine[igrid] ) continue;
 386     TList* padList = canKine[igrid]->GetListOfPrimitives();
 387     currName = canKine[igrid]->GetName();
 388     currName.Append(Form("_%s", basePlotName.Data()));
 389     can = new TCanvas(currName.Data(),currName.Data(),canKine[igrid]->GetWindowTopX(),igroup2*yshift,600,600);
 390     can->Divide(2,2);
 391     for ( Int_t ipad=0; ipad<padList->GetEntries(); ipad++ ) {
 392       TPad* pad = dynamic_cast<TPad*> (padList->At(ipad));
 393       if ( ! pad ) continue;
 394       TList* histoList = pad->GetListOfPrimitives();
 395       can->cd(ipad+1);
 396       for ( Int_t iobj=0; iobj<histoList->GetEntries(); iobj++ ) {
 397         currName = histoList->At(iobj)->GetName();
 398         if ( ! histoList->At(iobj)->InheritsFrom(TH1::Class()) || ! currName.Contains(fChargeKeys->At(1)->GetName()) ) continue;
 399         histoName = currName;
 400         histoName.ReplaceAll(fChargeKeys->At(1)->GetName(),"");
 401         histoName.Append(Form("_%s", basePlotName.Data()));
 402         currName.ReplaceAll(fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName());
 403         TH1* auxHisto = dynamic_cast<TH1*> (histoList->FindObject(currName.Data()));
 404         if ( ! auxHisto ) continue;
 405         TH1* histo = static_cast<TH1*> (histoList->At(iobj)->Clone(histoName.Data()));
 406         if ( plotChargeAsymmetry ) {
 407           histo->Add(auxHisto, -1.);
 408           // h2 + h1 = 2xh2 + (h1-h2)
 409           auxHisto->Add(auxHisto, histo, 2.);
 410         }
 411         histo->Divide(auxHisto);
 412         histo->SetMarkerStyle(20);
 413         TString axisTitle = plotChargeAsymmetry ? Form("(%s - %s) / (%s + %s)", fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName(), fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName()) : Form("%s / %s", fChargeKeys->At(1)->GetName(), fChargeKeys->At(0)->GetName());
 414         axisTitle.ReplaceAll("MuPlus","#mu^{+}");
 415         axisTitle.ReplaceAll("MuMinus","#mu^{-}");
 416         histo->GetYaxis()->SetTitle(axisTitle.Data());
 417         histo->SetStats(kFALSE);
 418         drawOpt = ( gPad->GetListOfPrimitives()->GetEntries() == 0 ) ? "e" : "esames";
 419         histo->Draw(drawOpt.Data());
 420       } // loop on histos
 421       gPad->Update();
 422     } // loop on pads
 423   } // loop on container steps
 424
 425
 426   //////////////////////
 427   // Event statistics //
 428   //////////////////////
 429   printf("\nTotal analyzed events:\n");
 430   TString evtSel = Form("trigger:%s", trigClassName.Data());
 431   fEventCounters->PrintSum(evtSel.Data());
 432   printf("Physics selected analyzed events:\n");
 433   evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
 434   fEventCounters->PrintSum(evtSel.Data());
 435
 436   TString countPhysSel = "any";
 437   if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
 438   else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
 439   countPhysSel.Prepend("selected:");
 440   printf("Analyzed events vs. centrality:\n");
 441   evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
 442   fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
 443
 444
 445   ///////////////////
 446   // Vertex method //
 447   ///////////////////
 448   if ( ! furtherOpt.Contains("VERTEX") ) return;
 449   Int_t firstMother = ( isMC ) ? 0 : kUnidentified;
 450   Int_t lastMother = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
 451   igroup1++;
 452   TH1* eventVertex = (TH1*)GetSum(physSel, minBiasTrig, centralityRange, "hIpVtx");
 453   if ( ! eventVertex ) return;
 454   Double_t minZ = -9.99, maxZ = 9.99;
 455   Double_t meanZ = 0., sigmaZ = 4.;
 456   Double_t nSigma = 2.;
 457   TString fitOpt = "R0S";
 458   Bool_t fixFitRange = kFALSE;
 459   TString fitFormula = Form("[0]+[1]*(x+[2])");
 460
 461   // Get vertex shape
 462   if ( eventVertex->GetSumw2N() == 0 ) eventVertex->Sumw2();
 463   Double_t eventVtxIntegral = eventVertex->Integral(0,eventVertex->GetNbinsX()+1); // Include under/overflow
 464   printf("Event vertex integral %.0f\n\n", eventVtxIntegral);
 465   if ( eventVtxIntegral <= 0. ) return;
 466   eventVertex->Scale(1./eventVtxIntegral);
 467   printf("\nFit MB vertex\n");
 468   eventVertex->Fit("gaus",fitOpt.Data(),"",minZ,maxZ);
 469   TF1* vtxFit = (TF1*)eventVertex->GetListOfFunctions()->FindObject("gaus");
 470   currName = "vtxIntegrated";
 471   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 472   can->SetLogy();
 473   eventVertex->Draw();
 474   vtxFit->Draw("same");
 475
 476
 477   enum {kRecoHF, kRecoBkg, kInputHF, kInputDecay, kRecoAll, kNrecoHistos};
 478   TString baseRecoName[kNrecoHistos] = {"RecoHF", "RecoBkg", "InputHF", "InputDecay", "RecoAll"};
 479   TArrayI sumMothers[kNrecoHistos];
 480   sumMothers[kRecoHF].Set(0);
 481   sumMothers[kRecoBkg].Set(0);
 482   sumMothers[kInputHF].Set(3);
 483   sumMothers[kInputHF][0] = kCharmMu;
 484   sumMothers[kInputHF][1] = kBeautyMu;
 485   sumMothers[kInputHF][2] = kQuarkoniumMu;
 486   sumMothers[kInputDecay].Set(1);
 487   sumMothers[kInputDecay][0] = kDecayMu;
 488   sumMothers[kRecoAll].Set(kNtrackSources);
 489   for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 490     sumMothers[kRecoAll][isrc] = isrc;
 491   }
 492
 493   meanZ = vtxFit->GetParameter(1);
 494   sigmaZ = vtxFit->GetParameter(2);
 495
 496   Double_t minZfit = ( fixFitRange ) ? minZ : meanZ - nSigma*sigmaZ;
 497   Double_t maxZfit = ( fixFitRange ) ? maxZ : meanZ + nSigma*sigmaZ;
 498
 499   TF1* fitFunc = new TF1("fitFunc", fitFormula.Data(), minZ, maxZ);
 500   fitFunc->SetLineColor(2);
 501   fitFunc->SetParNames("Line norm", "Line slope", "Free path");
 502   const Double_t kFreePath = 153.; // 150.; // 130.; // cm
 503   //fitFunc->SetParameters(0.,1.);
 504   fitFunc->FixParameter(2, kFreePath);
 505
 506   AliCFGridSparse* gridSparse = cfContainer->GetGrid(kStepReconstructed);
 507   TAxis* ptAxis = gridSparse->GetAxis(kHvarPt);
 508
 509   Double_t slope = 0.;
 510   Double_t limitNorm = 0., limitSlope = 0.;
 511   Int_t firstPtBin = 0, lastPtBin = 0;
 512
 513   gStyle->SetOptFit(1111);
 514
 515   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 516     igroup2++;
 517     SetSparseRange(gridSparse, kHvarThetaAbs, "", itheta+1, itheta+1, "USEBIN");
 518     SetSparseRange(gridSparse, kHvarPt, "", 1, ptAxis->GetNbins(), "USEBIN");
 519     TH1* recoHisto[kNrecoHistos];
 520     for ( Int_t ireco=0; ireco<kNrecoHistos; ++ireco ) {
 521       recoHisto[ireco] = gridSparse->Project(kHvarPt);
 522       histoName = Form("%sMuon_%s", baseRecoName[ireco].Data(), fThetaAbsKeys->At(itheta)->GetName());
 523       recoHisto[ireco]->SetName(histoName.Data());
 524       recoHisto[ireco]->SetTitle(histoName.Data());
 525       recoHisto[ireco]->Reset();
 526       recoHisto[ireco]->Sumw2();
 527       for ( Int_t isrc=0; isrc<sumMothers[ireco].GetSize(); ++isrc ) {
 528         SetSparseRange(gridSparse, kHvarMotherType, "", sumMothers[ireco][isrc]+1, sumMothers[ireco][isrc]+1, "USEBIN");
 529         TH1* auxHisto = gridSparse->Project(kHvarPt);
 530         recoHisto[ireco]->Add(auxHisto);
 531         delete auxHisto;
 532       }
 533     }
 534     SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
 535     Int_t currDraw = 0;
 536
 537     for ( Int_t ibinpt=0; ibinpt<=ptAxis->GetNbins(); ++ibinpt ) {
 538       firstPtBin = ibinpt;
 539       lastPtBin = ( ibinpt == 0 ) ? ptAxis->GetNbins() : ibinpt;
 540       SetSparseRange(gridSparse, kHvarPt, "", firstPtBin, lastPtBin, "USEBIN");
 541       TH1* histo = gridSparse->Project(kHvarVz);
 542       histo->SetName(Form("hVtx_%s_%s_ptBin%i", cfContainer->GetStepTitle(kStepReconstructed), fThetaAbsKeys->At(itheta)->GetName(), ibinpt));
 543       if ( histo->GetEntries() < 100. ) break;
 544       printf("\nFit %.2f < pt < %.2f (entries %g)\n", ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin), histo->GetEntries());
 545       histo->Divide(eventVertex);
 546       Double_t norm = histo->GetBinContent(histo->FindBin(0.));
 547       histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)");
 548       slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
 549                histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ );
 550
 551       if ( slope < 0. ) slope = norm/kFreePath;
 552
 553       // Try to fit twice: it fit fails the first time
 554       // set some limits on parameters
 555       for ( Int_t itry=0; itry<2; itry++ ) {
 556         fitFunc->SetParameter(0, norm);
 557         fitFunc->SetParameter(1, slope);
 558         if ( itry > 0 ) {
 559           limitNorm = 2.*histo->Integral();
 560           limitSlope = 2.*histo->Integral()/kFreePath;
 561           //fitFunc->SetParLimits(0, 0., limitNorm); // REMEMBER TO CHECK
 562           fitFunc->SetParLimits(1, 0., limitSlope); // REMEMBER TO CHECK
 563           printf("Norm 0. < %f < %f  slope  0. < %f < %f\n", norm, limitNorm, slope, limitSlope);
 564         }
 565         TFitResultPtr fitRes = histo->Fit(fitFunc, fitOpt.Data(), "", minZfit, maxZfit);
 566
 567 //      if ( gMinuit->fCstatu.Contains("CONVERGED") &&
 568         if ( ((Int_t)fitRes) == 0 &&
 569             fitFunc->GetParameter(0) > 0. &&
 570             fitFunc->GetParameter(1) > 0. )
 571           break;
 572         else if ( furtherOpt.Contains("REFIT") ) printf("Re-fit with limits\n");
 573         else {
 574           printf("Warning: fit problems !!!\n");
 575           break;
 576         }
 577       }
 578
 579       Double_t p0 = fitFunc->GetParameter(0);
 580       Double_t p0err = fitFunc->GetParError(0);
 581       Double_t p1 = fitFunc->GetParameter(1);
 582       Double_t p1err = fitFunc->GetParError(1);
 583
 584       Double_t nullVz = ( p1 != 0. ) ? -p0/p1 : 0.;
 585       Double_t nullVzErr = ( p0 != 0. && p1 != 0. ) ? TMath::Abs(nullVz) * TMath::Sqrt(p0err*p0err/(p0*p0) + p1err*p1err/(p1*p1) ) : 0.;
 586
 587       printf("Null value at %f +- %f\n", nullVz - kFreePath, nullVzErr);
 588
 589       recoHisto[kRecoHF]->SetBinContent(ibinpt, p0);
 590       recoHisto[kRecoHF]->SetBinError(ibinpt, p0err);
 591       recoHisto[kRecoBkg]->SetBinContent(ibinpt, ( kFreePath + meanZ ) * p1);
 592       recoHisto[kRecoBkg]->SetBinError(ibinpt, ( kFreePath + meanZ ) * p1err);
 593       if ( currDraw%4 == 0 ){
 594         currName = Form("vtx_%s_PtBin%i",fThetaAbsKeys->At(itheta)->GetName(), ibinpt);
 595         can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 596         can->Divide(2,2);
 597       }
 598       can->cd( currDraw%4 + 1 );
 599       can->SetLogy();
 600       histo->Draw();
 601       fitFunc->DrawCopy("same");
 602       currDraw++;
 603     } // loop on pt bins
 604     SetSparseRange(gridSparse, kHvarMotherType, "", firstMother+1, lastMother+1, "USEBIN");
 605     currName = Form("recoPt_%s",fThetaAbsKeys->At(itheta)->GetName());
 606     can = new TCanvas(currName.Data(),currName.Data(),(igroup1+1)*xshift,igroup2*yshift,600,600);
 607     TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8);
 608     drawOpt = "e";
 609     for ( Int_t ireco=0; ireco<kNrecoHistos-1; ++ireco ) {
 610       if ( recoHisto[ireco]->GetEntries() == 0. ) continue;
 611       TH1* ratio = (TH1*)recoHisto[ireco]->Clone(Form("%s_ratio", recoHisto[ireco]->GetName()));
 612       ratio->Divide(recoHisto[kRecoAll]);
 613       ratio->SetLineColor(srcColors[ireco]);
 614       ratio->SetMarkerColor(srcColors[ireco]);
 615       ratio->SetMarkerStyle(20+ireco);
 616       ratio->GetYaxis()->SetTitle("fraction of total");
 617       ratio->Draw(drawOpt.Data());
 618       leg->AddEntry(ratio,baseRecoName[ireco].Data(), "lp");
 619       drawOpt = "esame";
 620     }
 621     leg->Draw("same");
 622   } // loop on theta abs
 623 }