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