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