PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 AliAnalysisTaskSEDvsMultiplicity
  20 // AliAnalysisTaskSE for the D meson vs. multiplcity analysis
  21 // Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TClonesArray.h>
  25 #include <TCanvas.h>
  26 #include <TList.h>
  27 #include <TString.h>
  28 #include <TDatabasePDG.h>
  29 #include <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TH3F.h>
  32 #include <THnSparse.h>
  33 #include <TProfile.h>
  34 #include "AliAnalysisManager.h"
  35 #include "AliRDHFCuts.h"
  36 #include "AliRDHFCutsDplustoKpipi.h"
  37 #include "AliRDHFCutsDStartoKpipi.h"
  38 #include "AliRDHFCutsD0toKpi.h"
  39 #include "AliAODHandler.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODVertex.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODRecoDecayHF.h"
  44 #include "AliAODRecoCascadeHF.h"
  45 #include "AliAnalysisVertexingHF.h"
  46 #include "AliAnalysisTaskSE.h"
  47 #include "AliAnalysisTaskSEDvsMultiplicity.h"
  48 #include "AliNormalizationCounter.h"
  49 #include "AliVertexingHFUtils.h"
  50 #include "AliAODVZERO.h"
  51 #include "AliESDUtils.h"
  52 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
  53
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
  57 AliAnalysisTaskSE(),
  58   fOutput(0),
  59   fListCuts(0),
  60   fOutputCounters(0),
  61   fListProfiles(0),
  62   fHistNEvents(0),
  63   fHistNtrEta16vsNtrEta1EvSel(0),
  64   fHistNtrEta05vsNtrEta1EvSel(0),
  65   fHistNtrEta03vsNtrEta1EvSel(0),
  66   fHistNtrEtaV0AvsNtrEta1EvSel(0),
  67   fHistNtrEtaV0MvsNtrEta1EvSel(0),
  68   fHistNtrEtaV0AvsV0AEqEvSel(0),
  69   fHistNtrEtaV0MvsV0MEqEvSel(0),
  70   fHistNtrCorrEta1vsNtrRawEta1EvSel(0),
  71   fHistMultCorrvsMultRawEvSel(0),
  72   fHistNtrEta16vsNtrEta1EvWithCand(0),
  73   fHistNtrEta05vsNtrEta1EvWithCand(0),
  74   fHistNtrEta03vsNtrEta1EvWithCand(0),
  75   fHistNtrEtaV0AvsNtrEta1EvWithCand(0),
  76   fHistNtrEtaV0MvsNtrEta1EvWithCand(0),
  77   fHistNtrEtaV0AvsV0AEqEvWithCand(0),
  78   fHistNtrEtaV0MvsV0MEqEvWithCand(0),
  79   fHistNtrCorrEta1vsNtrRawEta1EvWithCand(0),
  80   fHistMultCorrvsMultRawEvWithCand(0),
  81   fHistNtrEta16vsNtrEta1EvWithD(0),
  82   fHistNtrEta05vsNtrEta1EvWithD(0),
  83   fHistNtrEta03vsNtrEta1EvWithD(0),
  84   fHistNtrEtaV0AvsNtrEta1EvWithD(0),
  85   fHistNtrEtaV0MvsNtrEta1EvWithD(0),
  86   fHistNtrEtaV0AvsV0AEqEvWithD(0),
  87   fHistNtrEtaV0MvsV0MEqEvWithD(0),
  88   fHistNtrCorrEta1vsNtrRawEta1EvWithD(0),
  89   fHistMultCorrvsMultRawEvWithD(0),
  90   fHistNtrVsZvtx(0),
  91   fHistNtrCorrVsZvtx(0),
  92   fHistNtrVsNchMC(0),
  93   fHistNtrCorrVsNchMC(0),
  94   fHistNtrVsNchMCPrimary(0),
  95   fHistNtrCorrVsNchMCPrimary(0),
  96   fHistNtrVsNchMCPhysicalPrimary(0),
  97   fHistNtrCorrVsNchMCPhysicalPrimary(0),
  98   fHistGenPrimaryParticlesInelGt0(0),
  99   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
 100   fHistNtrUnCorrPSSel(0),
 101   fHistNtrUnCorrPSTrigSel(0),
 102   fHistNtrUnCorrPSTrigPileUpSel(0),
 103   fHistNtrUnCorrPSTrigPileUpVtxSel(0),
 104   fHistNtrUnCorrPSTrigPileUpVtxContSel(0),
 105   fHistNtrUnCorrPSTrigPileUpVtxRangeSel(0),
 106   fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel(0),
 107   fHistNtrUnCorrEvSel(0),
 108   fHistNtrUnCorrEvWithCand(0),
 109   fHistNtrUnCorrEvWithD(0),
 110   fHistNtrCorrPSSel(0),
 111   fHistNtrCorrEvSel(0),
 112   fHistNtrCorrEvWithCand(0),
 113   fHistNtrCorrEvWithD(0),
 114   fPtVsMassVsMult(0),
 115   fPtVsMassVsMultNoPid(0),
 116   fPtVsMassVsMultUncorr(0),
 117   fPtVsMassVsMultPart(0),
 118   fPtVsMassVsMultAntiPart(0),
 119   fPtVsMassVsMultMC(0),
 120   fUpmasslimit(1.965),
 121   fLowmasslimit(1.765),
 122   fNMassBins(200),
 123   fRDCutsAnalysis(0),
 124   fCounter(0),
 125   fCounterU(0),
 126   fDoImpPar(kFALSE),
 127   fNImpParBins(400),
 128   fLowerImpPar(-2000.),
 129   fHigherImpPar(2000.),
 130   fReadMC(kFALSE),
 131   fMCOption(0),
 132   fisPPbData(kFALSE),
 133   fUseBit(kTRUE),
 134   fSubtractTrackletsFromDau(kFALSE),
 135   fKeepCorrPlots(kFALSE),
 136   fUseNchWeight(kFALSE),
 137   fHistoMCNch(0),
 138   fHistoMeasNch(0),
 139   fRefMult(9.26),
 140   fPdgMeson(411),
 141   fMultiplicityEstimator(kNtrk10),
 142   fMCPrimariesEstimator(kEta10),
 143   fDoVZER0ParamVertexCorr(0)
 144 {
 145    // Default constructor
 146   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 147   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 148 }
 149
 150 //________________________________________________________________________
 151 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts, Bool_t switchPPb):
 152   AliAnalysisTaskSE(name),
 153   fOutput(0),
 154   fListCuts(0),
 155   fOutputCounters(0),
 156   fListProfiles(0),
 157   fHistNEvents(0),
 158   fHistNtrEta16vsNtrEta1EvSel(0),
 159   fHistNtrEta05vsNtrEta1EvSel(0),
 160   fHistNtrEta03vsNtrEta1EvSel(0),
 161   fHistNtrEtaV0AvsNtrEta1EvSel(0),
 162   fHistNtrEtaV0MvsNtrEta1EvSel(0),
 163   fHistNtrEtaV0AvsV0AEqEvSel(0),
 164   fHistNtrEtaV0MvsV0MEqEvSel(0),
 165   fHistNtrCorrEta1vsNtrRawEta1EvSel(0),
 166   fHistMultCorrvsMultRawEvSel(0),
 167   fHistNtrEta16vsNtrEta1EvWithCand(0),
 168   fHistNtrEta05vsNtrEta1EvWithCand(0),
 169   fHistNtrEta03vsNtrEta1EvWithCand(0),
 170   fHistNtrEtaV0AvsNtrEta1EvWithCand(0),
 171   fHistNtrEtaV0MvsNtrEta1EvWithCand(0),
 172   fHistNtrEtaV0AvsV0AEqEvWithCand(0),
 173   fHistNtrEtaV0MvsV0MEqEvWithCand(0),
 174   fHistNtrCorrEta1vsNtrRawEta1EvWithCand(0),
 175   fHistMultCorrvsMultRawEvWithCand(0),
 176   fHistNtrEta16vsNtrEta1EvWithD(0),
 177   fHistNtrEta05vsNtrEta1EvWithD(0),
 178   fHistNtrEta03vsNtrEta1EvWithD(0),
 179   fHistNtrEtaV0AvsNtrEta1EvWithD(0),
 180   fHistNtrEtaV0MvsNtrEta1EvWithD(0),
 181   fHistNtrEtaV0AvsV0AEqEvWithD(0),
 182   fHistNtrEtaV0MvsV0MEqEvWithD(0),
 183   fHistNtrCorrEta1vsNtrRawEta1EvWithD(0),
 184   fHistMultCorrvsMultRawEvWithD(0),
 185   fHistNtrVsZvtx(0),
 186   fHistNtrCorrVsZvtx(0),
 187   fHistNtrVsNchMC(0),
 188   fHistNtrCorrVsNchMC(0),
 189   fHistNtrVsNchMCPrimary(0),
 190   fHistNtrCorrVsNchMCPrimary(0),
 191   fHistNtrVsNchMCPhysicalPrimary(0),
 192   fHistNtrCorrVsNchMCPhysicalPrimary(0),
 193   fHistGenPrimaryParticlesInelGt0(0),
 194   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
 195   fHistNtrUnCorrPSSel(0),
 196   fHistNtrUnCorrPSTrigSel(0),
 197   fHistNtrUnCorrPSTrigPileUpSel(0),
 198   fHistNtrUnCorrPSTrigPileUpVtxSel(0),
 199   fHistNtrUnCorrPSTrigPileUpVtxContSel(0),
 200   fHistNtrUnCorrPSTrigPileUpVtxRangeSel(0),
 201   fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel(0),
 202   fHistNtrUnCorrEvSel(0),
 203   fHistNtrUnCorrEvWithCand(0),
 204   fHistNtrUnCorrEvWithD(0),
 205   fHistNtrCorrPSSel(0),
 206   fHistNtrCorrEvSel(0),
 207   fHistNtrCorrEvWithCand(0),
 208   fHistNtrCorrEvWithD(0),
 209   fPtVsMassVsMult(0),
 210   fPtVsMassVsMultNoPid(0),
 211   fPtVsMassVsMultUncorr(0),
 212   fPtVsMassVsMultPart(0),
 213   fPtVsMassVsMultAntiPart(0),
 214   fPtVsMassVsMultMC(0),
 215   fUpmasslimit(1.965),
 216   fLowmasslimit(1.765),
 217   fNMassBins(200),
 218   fRDCutsAnalysis(cuts),
 219   fCounter(0),
 220   fCounterU(0),
 221   fDoImpPar(kFALSE),
 222   fNImpParBins(400),
 223   fLowerImpPar(-2000.),
 224   fHigherImpPar(2000.),
 225   fReadMC(kFALSE),
 226   fMCOption(0),
 227   fisPPbData(switchPPb),
 228   fUseBit(kTRUE),
 229   fSubtractTrackletsFromDau(kFALSE),
 230   fKeepCorrPlots(kFALSE),
 231   fUseNchWeight(kFALSE),
 232   fHistoMCNch(0),
 233   fHistoMeasNch(0),
 234   fRefMult(9.26),
 235   fPdgMeson(pdgMeson),
 236   fMultiplicityEstimator(kNtrk10),
 237   fMCPrimariesEstimator(kEta10),
 238   fDoVZER0ParamVertexCorr(0)
 239 {
 240   //
 241   // Standard constructor
 242   //
 243
 244   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 245   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 246   if(fPdgMeson==413){
 247     fNMassBins=200;
 248     SetMassLimits(0.12,0.2);
 249   }else{
 250     fNMassBins=200;
 251     SetMassLimits(fPdgMeson,0.1);
 252   }
 253   // Default constructor
 254    // Otput slot #1 writes into a TList container
 255   DefineOutput(1,TList::Class());  //My private output
 256   // Output slot #2 writes cut to private output
 257   DefineOutput(2,TList::Class());
 258   // Output slot #3 writes cut to private output
 259   DefineOutput(3,TList::Class());
 260   // Output slot #4 writes cut to private output
 261   DefineOutput(4,TList::Class());
 262 }
 263 //________________________________________________________________________
 264 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
 265 {
 266   //
 267   // Destructor
 268   //
 269   delete fOutput;
 270   delete fHistNEvents;
 271   delete fListCuts;
 272   delete fListProfiles;
 273   delete fRDCutsAnalysis;
 274   delete fCounter;
 275   delete fCounterU;
 276   for(Int_t i=0; i<4; i++) {
 277       if (fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i];
 278   }
 279
 280   for(Int_t i=0; i<5; i++){
 281     delete fHistMassPtImpPar[i];
 282   }
 283   if(fHistoMCNch) delete fHistoMCNch;
 284   if(fHistoMeasNch) delete fHistoMeasNch;
 285 }
 286
 287 //_________________________________________________________________
 288 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 289   // set invariant mass limits
 290   if(uplimit>lowlimit){
 291     fLowmasslimit = lowlimit;
 292     fUpmasslimit = uplimit;
 293   }else{
 294     AliError("Wrong mass limits: upper value should be larger than lower one");
 295   }
 296 }
 297 //_________________________________________________________________
 298 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
 299   // set invariant mass limits
 300   Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 301   SetMassLimits(mass-range,mass+range);
 302 }
 303 //________________________________________________________________________
 304 void AliAnalysisTaskSEDvsMultiplicity::Init(){
 305   //
 306   // Initialization
 307   //
 308   printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
 309
 310   if(fUseNchWeight && !fReadMC){ AliFatal("Nch weights can only be used in MC mode"); return; }
 311   if(fUseNchWeight && !fHistoMCNch){ AliFatal("Nch weights can only be used without histogram"); return; }
 312
 313   fListCuts=new TList();
 314   fListCuts->SetOwner();
 315   fListCuts->SetName("CutsList");
 316
 317
 318   if(fPdgMeson==411){
 319      AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
 320     copycut->SetName("AnalysisCutsDplus");
 321     fListCuts->Add(copycut);
 322   }else if(fPdgMeson==421){
 323     AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
 324     copycut->SetName("AnalysisCutsDzero");
 325     fListCuts->Add(copycut);
 326   }else if(fPdgMeson==413){
 327      AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
 328     copycut->SetName("AnalysisCutsDStar");
 329     fListCuts->Add(copycut);
 330   }
 331   PostData(2,fListCuts);
 332
 333   fListProfiles = new TList();
 334   fListProfiles->SetOwner();
 335   TString period[4];
 336   Int_t nProfiles=4;
 337   if (fisPPbData) {period[0]="LHC13b"; period[1]="LHC13c"; nProfiles = 2;}
 338   else {period[0]="LHC10b"; period[1]="LHC10c"; period[2]="LHC10d"; period[3]="LHC10e"; nProfiles = 4;}
 339
 340   for(Int_t i=0; i<nProfiles; i++){
 341     if(fMultEstimatorAvg[i]){
 342       TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
 343       hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
 344       fListProfiles->Add(hprof);
 345     }
 346   }
 347   PostData(4,fListProfiles);
 348
 349   return;
 350 }
 351
 352 //________________________________________________________________________
 353 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
 354 {
 355   // Create the output container
 356   //
 357   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
 358
 359   // Several histograms are more conveniently managed in a TList
 360   fOutput = new TList();
 361   fOutput->SetOwner();
 362   fOutput->SetName("OutputHistos");
 363
 364   Int_t nMultBins = 200;
 365   Float_t firstMultBin = -0.5;
 366   Float_t lastMultBin = 199.5;
 367   Int_t nMultBinsNtrk = nMultBins;
 368   Float_t lastMultBinNtrk = lastMultBin;
 369   Int_t nMultBinsV0 = 400;
 370   Float_t lastMultBinV0 = 799.5;
 371   const char *estimatorName="tracklets";
 372   if(fisPPbData) {
 373     nMultBinsNtrk = 375;
 374     lastMultBinNtrk = 374.5;
 375     nMultBins = nMultBinsNtrk;
 376     lastMultBin = lastMultBinNtrk;
 377   }
 378   if(fMultiplicityEstimator==kVZERO || fMultiplicityEstimator==kVZEROA || fMultiplicityEstimator==kVZEROEq || fMultiplicityEstimator==kVZEROAEq) {
 379     nMultBins = nMultBinsV0;
 380     lastMultBin = lastMultBinV0;
 381     estimatorName = "vzero";
 382   }
 383
 384   fHistNtrUnCorrPSSel = new TH1F("hNtrUnCorrPSSel",Form("Uncorrected %s multiplicity for PS selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 385   fHistNtrUnCorrPSTrigSel = new TH1F("hNtrUnCorrPSTrigSel",Form("Uncorrected %s multiplicity for PS + trigger name selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 386   fHistNtrUnCorrPSTrigPileUpSel = new TH1F("hNtrUnCorrPSTrigPileUpSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 387   fHistNtrUnCorrPSTrigPileUpVtxSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 388   fHistNtrUnCorrPSTrigPileUpVtxContSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxContSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex-contrib selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 389   fHistNtrUnCorrPSTrigPileUpVtxRangeSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxRangeSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex-contrib-range selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 390   fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel = new TH1F("hNtrUnCorrPSTrigPileUpVtxRangeCentrSel",Form("Uncorrected %s multiplicity for PS + trigger name + pileup + with-vertex-contrib-range + centrality selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 391   fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel",Form("Uncorrected %s multiplicity for selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 392   fHistNtrUnCorrEvWithCand = new TH1F("hNtrUnCorrEvWithCand",Form("Uncorrected %s multiplicity for events with D candidates; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);// Total multiplicity
 393   fHistNtrUnCorrEvWithD = new TH1F("hNtrUnCorrEvWithD",Form("Uncorrected %s multiplicity for events with D in mass region ; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin); //
 394   fHistNtrCorrPSSel = new TH1F("hNtrCorrPSSel",Form("Corrected %s multiplicity for PS selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 395   fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel",Form("Corrected %s multiplicity for selected events; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);
 396   fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", Form("%s multiplicity for events with D candidates; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin);// Total multiplicity
 397   fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", Form("%s multiplicity for events with D in mass region ; %s ; Entries",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin); //
 398
 399   if(fKeepCorrPlots){
 400     fHistNtrEta16vsNtrEta1EvSel = new TH2F("hNtrEta16vsNtrEta1EvSel","Uncorrected Eta1.6 vs Eta1.0 (events selected); Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
 401     fHistNtrEta05vsNtrEta1EvSel = new TH2F("hNtrEta05vsNtrEta1EvSel","Uncorrected Eta0.5 vs Eta1.0 (events selected); Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.5 vs eta 1.0 histogram
 402     fHistNtrEta03vsNtrEta1EvSel = new TH2F("hNtrEta03vsNtrEta1EvSel","Uncorrected Eta0.3 vs Eta1.0 (events selected); Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.3 vs eta 1.0 histogram
 403     fHistNtrEtaV0AvsNtrEta1EvSel = new TH2F("hNtrEtaV0AvsNtrEta1EvSel","Uncorrected Eta-V0A vs Eta1.0 (events selected); Ntracklets #eta<1.0; Multiplicity V0A",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0A vs eta 1.0 histogram
 404     fHistNtrEtaV0MvsNtrEta1EvSel = new TH2F("hNtrEtaV0MvsNtrEta1EvSel","Uncorrected Eta-V0M vs Eta1.0 (events selected); Ntracklets #eta<1.0; Multiplicity V0A+V0C",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0M vs eta 1.0 histogram
 405     fHistNtrEtaV0AvsV0AEqEvSel = new TH2F("hNtrEtaV0AvsV0AEqEvSel","Corrected V0A vs corrected V0A-Equalized (events selected); Vzero-A; Vzero-A Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0A - V0Aeq
 406     fHistNtrEtaV0MvsV0MEqEvSel = new TH2F("hNtrEtaV0MvsV0MEqEvSel","Corrected V0M vs corrected V0M-Equalized (events selected); Vzero-M; Vzero-M Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0M - V0Meq
 407     fHistNtrCorrEta1vsNtrRawEta1EvSel = new TH2F("hNtrCorrEta1vsNtrRawEta1EvSel","Corrected Eta1 vs Eta1.0 (events selected); Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
 408     fHistMultCorrvsMultRawEvSel = new TH2F("hMultCorrvsMultRawEvSel",Form("Corrected multiplicity vs uncorrected multiplicity (events selected); %s corrected; %s",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); // corrected vs uncorrected multiplicity
 409
 410     fHistNtrEta16vsNtrEta1EvWithCand = new TH2F("hNtrEta16vsNtrEta1EvWithCand","Uncorrected Eta1.6 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
 411     fHistNtrEta05vsNtrEta1EvWithCand = new TH2F("hNtrEta05vsNtrEta1EvWithCand","Uncorrected Eta0.5 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.5 vs eta 1.0 histogram
 412     fHistNtrEta03vsNtrEta1EvWithCand = new TH2F("hNtrEta03vsNtrEta1EvWithCand","Uncorrected Eta0.3 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.3 vs eta 1.0 histogram
 413     fHistNtrEtaV0AvsNtrEta1EvWithCand = new TH2F("hNtrEtaV0AvsNtrEta1EvWithCand","Uncorrected Eta-V0A vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Multiplicity V0A",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0A vs eta 1.0 histogram
 414     fHistNtrEtaV0MvsNtrEta1EvWithCand = new TH2F("hNtrEtaV0MvsNtrEta1EvWithCand","Uncorrected Eta-V0M vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0; Multiplicity V0A+V0C",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0M vs eta 1.0 histogram
 415     fHistNtrEtaV0AvsV0AEqEvWithCand = new TH2F("hNtrEtaV0AvsV0AEqEvWithCand","Corrected V0A vs corrected V0A-Equalized (events selected with a D candidate); Vzero-A; Vzero-A Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0A - V0Aeq
 416     fHistNtrEtaV0MvsV0MEqEvWithCand = new TH2F("hNtrEtaV0MvsV0MEqEvWithCand","Corrected V0M vs corrected V0M-Equalized (events selected with a D candidate); Vzero-M; Vzero-M Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0M - V0Meq
 417     fHistNtrCorrEta1vsNtrRawEta1EvWithCand = new TH2F("hNtrCorrEta1vsNtrRawEta1EvWithCand","Corrected Eta1 vs Eta1.0 (events selected with a D candidate); Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
 418     fHistMultCorrvsMultRawEvWithCand = new TH2F("hMultCorrvsMultRawEvWithCand",Form("Corrected multiplicity vs uncorrected multiplicity (events selected) with a D candidate; %s corrected; %s",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); // corrected vs uncorrected multiplicity
 419
 420
 421     fHistNtrEta16vsNtrEta1EvWithD = new TH2F("hNtrEta16vsNtrEta1EvWithD","Uncorrected Eta1.6 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
 422     fHistNtrEta05vsNtrEta1EvWithD = new TH2F("hNtrEta05vsNtrEta1EvWithD","Uncorrected Eta0.5 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Ntracklets #eta<0.5",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.5 vs eta 1.0 histogram
 423     fHistNtrEta03vsNtrEta1EvWithD = new TH2F("hNtrEta03vsNtrEta1EvWithD","Uncorrected Eta0.3 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Ntracklets #eta<0.3",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 0.3 vs eta 1.0 histogram
 424     fHistNtrEtaV0AvsNtrEta1EvWithD = new TH2F("hNtrEtaV0AvsNtrEta1EvWithD","Uncorrected Eta-V0A vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Multiplicity V0A",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0A vs eta 1.0 histogram
 425     fHistNtrEtaV0MvsNtrEta1EvWithD = new TH2F("hNtrEtaV0MvsNtrEta1EvWithD","Uncorrected Eta-V0M vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0; Multiplicity V0A+V0C",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsV0,firstMultBin,lastMultBinV0); //eta V0M vs eta 1.0 histogram
 426     fHistNtrEtaV0AvsV0AEqEvWithD = new TH2F("hNtrEtaV0AvsV0AEqEvWithD","Corrected V0A vs corrected V0A-Equalized (events selected with D in mass range); Vzero-A; Vzero-A Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0A - V0Aeq
 427     fHistNtrEtaV0MvsV0MEqEvWithD = new TH2F("hNtrEtaV0MvsV0MEqEvWithD","Corrected V0M vs corrected V0M-Equalized (events selected with D in mass range); Vzero-M; Vzero-M Equalized",nMultBinsV0,firstMultBin,lastMultBinV0,nMultBinsV0,firstMultBin,lastMultBinV0); // comparison V0M - V0Meq
 428     fHistNtrCorrEta1vsNtrRawEta1EvWithD = new TH2F("hNtrCorrEta1vsNtrRawEta1EvWithD","Corrected Eta1 vs Eta1.0 (events selected with D in mass range); Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBinsNtrk,firstMultBin,lastMultBinNtrk,nMultBinsNtrk,firstMultBin,lastMultBinNtrk); //eta 1.6 vs eta 1.0 histogram
 429     fHistMultCorrvsMultRawEvWithD = new TH2F("hMultCorrvsMultRawEvWithD",Form("Corrected multiplicity vs uncorrected multiplicity (events selected with D in mass range); %s corrected; %s",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); // corrected vs uncorrected multiplicity
 430
 431   }
 432   fHistNtrVsZvtx = new TH2F("hNtrVsZvtx",Form("N%s vs VtxZ; VtxZ;N_{%s};",estimatorName,estimatorName),300,-15,15,nMultBins,firstMultBin,lastMultBin); //
 433   fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx",Form("N%s vs VtxZ; VtxZ;N_{%s};",estimatorName,estimatorName),300,-15,15,nMultBins,firstMultBin,lastMultBin); //
 434
 435   fHistNtrVsNchMC = new TH2F("hNtrVsNchMC",Form("N%s vs NchMC; Nch;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 436   fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC",Form("N%s vs Nch; Nch;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 437
 438   fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary",Form("N%s vs Nch (Primary); Nch (Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 439   fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary",Form("N%s vs Nch (Primary); Nch(Primary) ;N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 440
 441   fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary",Form("N%s vs Nch (Physical Primary); Nch (Physical Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 442   fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary",Form("N%s vs Nch (Physical Primary); Nch (Physical Primary);N_{%s};",estimatorName,estimatorName),nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 443
 444   fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",nMultBins,firstMultBin,lastMultBin);
 445
 446   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin);
 447
 448   fHistNtrUnCorrPSSel->Sumw2();
 449   fHistNtrUnCorrPSTrigSel->Sumw2();
 450   fHistNtrUnCorrPSTrigPileUpSel->Sumw2();
 451   fHistNtrUnCorrPSTrigPileUpVtxSel->Sumw2();
 452   fHistNtrUnCorrPSTrigPileUpVtxContSel->Sumw2();
 453   fHistNtrUnCorrPSTrigPileUpVtxRangeSel->Sumw2();
 454   fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel->Sumw2();
 455   fHistNtrUnCorrEvSel->Sumw2();
 456   fHistNtrUnCorrEvWithCand->Sumw2();
 457   fHistNtrUnCorrEvWithD->Sumw2();
 458   fHistNtrCorrPSSel->Sumw2();
 459   fHistNtrCorrEvSel->Sumw2();
 460   fHistNtrCorrEvWithCand->Sumw2();
 461   fHistNtrCorrEvWithD->Sumw2();
 462   fHistGenPrimaryParticlesInelGt0->Sumw2();
 463   fOutput->Add(fHistNtrUnCorrPSSel);
 464   fOutput->Add(fHistNtrUnCorrPSTrigSel);
 465   fOutput->Add(fHistNtrUnCorrPSTrigPileUpSel);
 466   fOutput->Add(fHistNtrUnCorrPSTrigPileUpVtxSel);
 467   fOutput->Add(fHistNtrUnCorrPSTrigPileUpVtxContSel);
 468   fOutput->Add(fHistNtrUnCorrPSTrigPileUpVtxRangeSel);
 469   fOutput->Add(fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel);
 470   fOutput->Add(fHistNtrUnCorrEvSel);
 471   fOutput->Add(fHistNtrUnCorrEvWithCand);
 472   fOutput->Add(fHistNtrUnCorrEvWithD);
 473   fOutput->Add(fHistNtrCorrPSSel);
 474   fOutput->Add(fHistNtrCorrEvSel);
 475   fOutput->Add(fHistNtrCorrEvWithCand);
 476   fOutput->Add(fHistNtrCorrEvWithD);
 477   if(fKeepCorrPlots){
 478     fOutput->Add(fHistNtrEta16vsNtrEta1EvSel);
 479     fOutput->Add(fHistNtrEta05vsNtrEta1EvSel);
 480     fOutput->Add(fHistNtrEta03vsNtrEta1EvSel);
 481     fOutput->Add(fHistNtrEtaV0AvsNtrEta1EvSel);
 482     fOutput->Add(fHistNtrEtaV0MvsNtrEta1EvSel);
 483     fOutput->Add(fHistNtrEtaV0AvsV0AEqEvSel);
 484     fOutput->Add(fHistNtrEtaV0MvsV0MEqEvSel);
 485     fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1EvSel);
 486     fOutput->Add(fHistMultCorrvsMultRawEvSel);
 487     fOutput->Add(fHistNtrEta16vsNtrEta1EvWithCand);
 488     fOutput->Add(fHistNtrEta05vsNtrEta1EvWithCand);
 489     fOutput->Add(fHistNtrEta03vsNtrEta1EvWithCand);
 490     fOutput->Add(fHistNtrEtaV0AvsNtrEta1EvWithCand);
 491     fOutput->Add(fHistNtrEtaV0MvsNtrEta1EvWithCand);
 492     fOutput->Add(fHistNtrEtaV0AvsV0AEqEvWithCand);
 493     fOutput->Add(fHistNtrEtaV0MvsV0MEqEvWithCand);
 494     fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1EvWithCand);
 495     fOutput->Add(fHistMultCorrvsMultRawEvWithCand);
 496     fOutput->Add(fHistNtrEta16vsNtrEta1EvWithD);
 497     fOutput->Add(fHistNtrEta05vsNtrEta1EvWithD);
 498     fOutput->Add(fHistNtrEta03vsNtrEta1EvWithD);
 499     fOutput->Add(fHistNtrEtaV0AvsNtrEta1EvWithD);
 500     fOutput->Add(fHistNtrEtaV0MvsNtrEta1EvWithD);
 501     fOutput->Add(fHistNtrEtaV0AvsV0AEqEvWithD);
 502     fOutput->Add(fHistNtrEtaV0MvsV0MEqEvWithD);
 503     fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1EvWithD);
 504     fOutput->Add(fHistMultCorrvsMultRawEvWithD);
 505   }
 506   fOutput->Add(fHistNtrVsZvtx);
 507   fOutput->Add(fHistNtrCorrVsZvtx);
 508
 509   fOutput->Add(fHistNtrVsNchMC);
 510   fOutput->Add(fHistNtrCorrVsNchMC);
 511   fOutput->Add(fHistNtrVsNchMCPrimary);
 512   fOutput->Add(fHistNtrCorrVsNchMCPrimary);
 513   fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
 514   fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
 515   fOutput->Add(fHistGenPrimaryParticlesInelGt0);
 516   fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
 517
 518
 519   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
 520   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
 521   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
 522   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
 523   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
 524   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
 525   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
 526   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
 527   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
 528   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 529   fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
 530   fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
 531   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 532   fHistNEvents->Sumw2();
 533   fHistNEvents->SetMinimum(0);
 534   fOutput->Add(fHistNEvents);
 535
 536   fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 537
 538   fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 539
 540   fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 541
 542   fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 543
 544   fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 545
 546   fPtVsMassVsMultMC=new TH3F("hPtVsMassvsMultMC", "D true candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 547
 548   fOutput->Add(fPtVsMassVsMult);
 549   fOutput->Add(fPtVsMassVsMultUncorr);
 550   fOutput->Add(fPtVsMassVsMultNoPid);
 551   fOutput->Add(fPtVsMassVsMultPart);
 552   fOutput->Add(fPtVsMassVsMultAntiPart);
 553   fOutput->Add(fPtVsMassVsMultMC);
 554
 555   if(fDoImpPar) CreateImpactParameterHistos();
 556
 557   fCounter = new AliNormalizationCounter("NormCounterCorrMult");
 558   fCounter->SetStudyMultiplicity(kTRUE,1.);
 559   fCounter->Init();
 560
 561   fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
 562   fCounterU->SetStudyMultiplicity(kTRUE,1.);
 563   fCounterU->Init();
 564
 565   fOutputCounters = new TList();
 566   fOutputCounters->SetOwner();
 567   fOutputCounters->SetName("OutputCounters");
 568   fOutputCounters->Add(fCounter);
 569   fOutputCounters->Add(fCounterU);
 570
 571   PostData(1,fOutput);
 572   PostData(2,fListCuts);
 573   PostData(3,fOutputCounters);
 574   PostData(4,fListProfiles);
 575
 576   if(fUseNchWeight) CreateMeasuredNchHisto();
 577
 578   return;
 579 }
 580
 581
 582
 583 //________________________________________________________________________
 584 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
 585 {
 586   // Execute analysis for current event:
 587   // heavy flavor candidates association to MC truth
 588
 589   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 590
 591   //  AliAODTracklets* tracklets = aod->GetTracklets();
 592   //Int_t ntracklets = tracklets->GetNumberOfTracklets();
 593
 594
 595   TClonesArray *arrayCand = 0;
 596   TString arrayName="";
 597   UInt_t pdgDau[3];
 598   Int_t nDau=0;
 599   Int_t selbit=0;
 600   if(fPdgMeson==411){
 601     arrayName="Charm3Prong";
 602     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
 603     nDau=3;
 604     selbit=AliRDHFCuts::kDplusCuts;
 605   }else if(fPdgMeson==421){
 606     arrayName="D0toKpi";
 607     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
 608     nDau=2;
 609     selbit=AliRDHFCuts::kD0toKpiCuts;
 610   }else if(fPdgMeson==413){
 611     arrayName="Dstar";
 612     pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=0; // Quoting here D0 daughters (D* ones on another variable later)
 613     nDau=2;
 614     selbit=AliRDHFCuts::kDstarCuts;
 615   }
 616
 617   if(!aod && AODEvent() && IsStandardAOD()) {
 618     // In case there is an AOD handler writing a standard AOD, use the AOD
 619     // event in memory rather than the input (ESD) event.
 620     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 621      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 622      // have to taken from the AOD event hold by the AliAODExtension
 623     AliAODHandler* aodHandler = (AliAODHandler*)
 624       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 625     if(aodHandler->GetExtensions()) {
 626       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 627       AliAODEvent *aodFromExt = ext->GetAOD();
 628       arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
 629     }
 630   } else if(aod) {
 631     arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
 632   }
 633
 634   if(!aod || !arrayCand) {
 635     printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
 636     return;
 637   }
 638
 639   if(fisPPbData && fReadMC){
 640     Int_t runnumber = aod->GetRunNumber();
 641     if(aod->GetTriggerMask()==0 &&
 642        (runnumber>=195344 && runnumber<=195677)){
 643       AliDebug(3,"Event rejected because of null trigger mask");
 644       return;
 645     }
 646   }
 647
 648
 649   // fix for temporary bug in ESDfilter
 650   // the AODs with null vertex pointer didn't pass the PhysSel
 651   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 652
 653   // Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 654   // Int_t countTreta03=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.3,0.3);
 655   // Int_t countTreta05=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-0.5,0.5);
 656   // Int_t countTreta16=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
 657   Int_t countTreta1=0, countTreta03=0, countTreta05=0, countTreta16=0;
 658   AliAODTracklets* tracklets=aod->GetTracklets();
 659   Int_t nTr=tracklets->GetNumberOfTracklets();
 660   for(Int_t iTr=0; iTr<nTr; iTr++){
 661     Double_t theta=tracklets->GetTheta(iTr);
 662     Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
 663     if(eta>-0.3 && eta<0.3) countTreta03++;
 664     if(eta>-0.5 && eta<0.5) countTreta05++;
 665     if(eta>-1.0 && eta<1.0) countTreta1++;
 666     if(eta>-1.6 && eta<1.6) countTreta16++;
 667   }
 668
 669
 670   Int_t vzeroMult=0, vzeroMultA=0, vzeroMultC=0;
 671   Int_t vzeroMultEq=0, vzeroMultAEq=0, vzeroMultCEq=0;
 672   AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
 673   if(vzeroAOD) {
 674     vzeroMultA = static_cast<Int_t>(vzeroAOD->GetMTotV0A());
 675     vzeroMultC = static_cast<Int_t>(vzeroAOD->GetMTotV0C());
 676     vzeroMult = vzeroMultA + vzeroMultC;
 677     vzeroMultAEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity(aod));
 678     vzeroMultCEq = static_cast<Int_t>(AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity(aod));
 679     vzeroMultEq = vzeroMultAEq + vzeroMultCEq;
 680   }
 681
 682   Int_t countMult = countTreta1;
 683   if(fMultiplicityEstimator==kNtrk03) { countMult = countTreta03; }
 684   else if(fMultiplicityEstimator==kNtrk05) { countMult = countTreta05; }
 685   else if(fMultiplicityEstimator==kNtrk10to16) { countMult = countTreta16 - countTreta1; }
 686   else if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
 687   else if(fMultiplicityEstimator==kVZEROA) { countMult = vzeroMultA; }
 688   else if(fMultiplicityEstimator==kVZEROEq) { countMult = vzeroMultEq; }
 689   else if(fMultiplicityEstimator==kVZEROAEq) { countMult = vzeroMultAEq; }
 690
 691
 692   fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countMult);
 693   fHistNEvents->Fill(0); // count event
 694
 695   Double_t countTreta1corr=countTreta1;
 696   Double_t countCorr=countMult;
 697   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 698   // In case of VZERO multiplicity, consider the zvtx correction flag
 699   //  fDoVZER0ParamVertexCorr: 0= none, 1= usual d2h, 2=AliESDUtils
 700   Bool_t isDataDrivenZvtxCorr=kTRUE;
 701   Bool_t isVtxOk=kFALSE;
 702   Int_t vzeroMultACorr=vzeroMultA, vzeroMultCCorr=vzeroMultC, vzeroMultCorr=vzeroMult;
 703   Int_t vzeroMultAEqCorr=vzeroMultAEq, vzeroMultCEqCorr=vzeroMultCEq, vzeroMultEqCorr=vzeroMultEq;
 704   if(vtx1){
 705     if(vtx1->GetNContributors()>0){
 706       fHistNEvents->Fill(1);
 707       isVtxOk=kTRUE;
 708     }
 709   }
 710   if(isVtxOk){
 711     if( (fMultiplicityEstimator==kVZERO) || (fMultiplicityEstimator==kVZEROA) ||
 712         (fMultiplicityEstimator==kVZEROEq) || (fMultiplicityEstimator==kVZEROAEq) ){
 713       if(fDoVZER0ParamVertexCorr==0){
 714         // do not correct
 715         isDataDrivenZvtxCorr=kFALSE;
 716       } else if (fDoVZER0ParamVertexCorr==2){
 717         // use AliESDUtils correction
 718         Float_t zvtx = vtx1->GetZ();
 719         isDataDrivenZvtxCorr=kFALSE;
 720         vzeroMultACorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultA,zvtx));
 721         vzeroMultCCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0C(vzeroMultC,zvtx));
 722         vzeroMultCorr = vzeroMultACorr + vzeroMultCCorr;
 723         vzeroMultAEqCorr = static_cast<Int_t>(AliESDUtils::GetCorrV0A(vzeroMultAEq,zvtx));
 724         vzeroMultCEqCorr =static_cast<Int_t>( AliESDUtils::GetCorrV0C(vzeroMultCEq,zvtx));
 725         vzeroMultEqCorr = vzeroMultAEqCorr + vzeroMultCEqCorr;
 726         if(fMultiplicityEstimator==kVZERO) { countCorr = vzeroMultCorr; }
 727         else if(fMultiplicityEstimator==kVZEROA) { countCorr = vzeroMultACorr; }
 728         else if(fMultiplicityEstimator==kVZEROEq) { countCorr = vzeroMultEqCorr; }
 729         else if(fMultiplicityEstimator==kVZEROAEq) { countCorr = vzeroMultAEqCorr; }
 730       }
 731     }
 732   }
 733   // Data driven multiplicity z-vertex correction
 734   if(isVtxOk && isDataDrivenZvtxCorr){
 735     TProfile* estimatorAvg = GetEstimatorHistogram(aod);
 736     if(estimatorAvg){
 737       countTreta1corr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult));
 738       // vzeroMultACorr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,vzeroMultA,vtx1->GetZ(),fRefMult));
 739       // vzeroMultCorr= vzeroMultACorr + static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,vzeroMultC,vtx1->GetZ(),fRefMult));
 740       // vzeroMultAEqCorr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,vzeroMultAEq,vtx1->GetZ(),fRefMult));
 741       // vzeroMultEqCorr= vzeroMultAEqCorr + static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,vzeroMultCEq,vtx1->GetZ(),fRefMult));
 742       countCorr=static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,vtx1->GetZ(),fRefMult));
 743     }
 744   }
 745
 746
 747   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 748
 749   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
 750   if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
 751   if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
 752   if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
 753
 754   Bool_t isEvPSRejected = fRDCutsAnalysis->IsEventRejectedDuePhysicsSelection();
 755   Bool_t isEvTrigNameRejected = fRDCutsAnalysis->IsEventRejectedDueToTrigger();
 756   Bool_t isEvPileUpRejected = fRDCutsAnalysis->IsEventRejectedDueToPileup();
 757   Bool_t isEvNoVtxRejected = fRDCutsAnalysis->IsEventRejectedDueToNotRecoVertex();
 758   Bool_t isEvVtxContribRejected = fRDCutsAnalysis->IsEventRejectedDueToVertexContributors();
 759   Bool_t isEvVtxRangeRejected= fRDCutsAnalysis->IsEventRejectedDueToZVertexOutsideFiducialRegion();
 760   Bool_t isEvCentralityRejected = fRDCutsAnalysis->IsEventRejectedDueToCentrality();
 761   if(!isEvPSRejected){
 762     fHistNtrUnCorrPSSel->Fill(countMult);
 763     fHistNtrCorrPSSel->Fill(countCorr);
 764     if(!isEvTrigNameRejected){
 765       fHistNtrUnCorrPSTrigSel->Fill(countMult);
 766       if(!isEvPileUpRejected){
 767         fHistNtrUnCorrPSTrigPileUpSel->Fill(countMult);
 768         if(!isEvNoVtxRejected){
 769           fHistNtrUnCorrPSTrigPileUpVtxSel->Fill(countMult);
 770           if(!isEvVtxContribRejected){
 771             fHistNtrUnCorrPSTrigPileUpVtxContSel->Fill(countMult);
 772             if(!isEvVtxRangeRejected){
 773               fHistNtrUnCorrPSTrigPileUpVtxRangeSel->Fill(countMult);
 774               if(!isEvCentralityRejected){
 775                 fHistNtrUnCorrPSTrigPileUpVtxRangeCentrSel->Fill(countMult);
 776               }
 777             }
 778           }
 779         }
 780       }
 781     }
 782   }
 783
 784   if(!isEvSel)return;
 785   if(fKeepCorrPlots){
 786     fHistNtrEta16vsNtrEta1EvSel->Fill(countTreta1,countTreta16);
 787     fHistNtrEta05vsNtrEta1EvSel->Fill(countTreta1,countTreta05);
 788     fHistNtrEta03vsNtrEta1EvSel->Fill(countTreta1,countTreta03);
 789     fHistNtrEtaV0AvsNtrEta1EvSel->Fill(countTreta1,vzeroMultA);
 790     fHistNtrEtaV0MvsNtrEta1EvSel->Fill(countTreta1,vzeroMult);
 791     fHistNtrEtaV0AvsV0AEqEvSel->Fill(vzeroMultA,vzeroMultAEq);
 792     fHistNtrEtaV0MvsV0MEqEvSel->Fill(vzeroMult,vzeroMultEq);
 793     fHistNtrCorrEta1vsNtrRawEta1EvSel->Fill(countTreta1,countTreta1corr);
 794     fHistMultCorrvsMultRawEvSel->Fill(countMult,countCorr);
 795   }
 796   if(vtx1){
 797     fHistNtrVsZvtx->Fill(vtx1->GetZ(),countMult);
 798     fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countCorr);
 799   }
 800
 801   TClonesArray *arrayMC=0;
 802   AliAODMCHeader *mcHeader=0;
 803
 804   Double_t nchWeight=1.0;
 805
 806   // load MC particles
 807   if(fReadMC){
 808
 809     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 810     if(!arrayMC) {
 811       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
 812       return;
 813     }
 814     // load MC header
 815     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 816     if(!mcHeader) {
 817       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
 818       return;
 819      }
 820
 821
 822     // Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
 823     // Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
 824     // Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
 825     //
 826     Int_t nChargedMCEta10=0, nChargedMCEta03=0, nChargedMCEta05=0, nChargedMCEta16=0, nChargedMCEtam37tm17=0, nChargedMCEta28t51=0;
 827     Int_t nChargedMCPrimaryEta10=0, nChargedMCPrimaryEta03=0, nChargedMCPrimaryEta05=0, nChargedMCPrimaryEta16=0, nChargedMCPrimaryEtam37tm17=0, nChargedMCPrimaryEta28t51=0;
 828     Int_t nChargedMCPhysicalPrimaryEta10=0, nChargedMCPhysicalPrimaryEta03=0, nChargedMCPhysicalPrimaryEta05=0, nChargedMCPhysicalPrimaryEta16=0, nChargedMCPhysicalPrimaryEtam37tm17=0, nChargedMCPhysicalPrimaryEta28t51=0;
 829     for(Int_t i=0; i<arrayMC->GetEntriesFast(); i++){
 830       AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
 831       Int_t charge = part->Charge();
 832       Double_t eta = part->Eta();
 833       Bool_t isPrim = part->IsPrimary();
 834       Bool_t isPhysPrim = part->IsPhysicalPrimary();
 835       if(charge!=0) {
 836         if(eta>-0.3 && eta< 0.3) {
 837           nChargedMCEta03++;
 838           if(isPrim) nChargedMCPrimaryEta03++;
 839           if(isPhysPrim) nChargedMCPhysicalPrimaryEta03++;
 840         }
 841         if(eta>-0.5 && eta< 0.5) {
 842           nChargedMCEta05++;
 843           if(isPrim) nChargedMCPrimaryEta05++;
 844           if(isPhysPrim) nChargedMCPhysicalPrimaryEta05++;
 845         }
 846         if(eta>-1.0 && eta< 1.0) {
 847           nChargedMCEta10++;
 848           if(isPrim) nChargedMCPrimaryEta10++;
 849           if(isPhysPrim) nChargedMCPhysicalPrimaryEta10++;
 850         }
 851         if(eta>-1.6 && eta< 1.6) {
 852           nChargedMCEta16++;
 853           if(isPrim) nChargedMCPrimaryEta16++;
 854           if(isPhysPrim) nChargedMCPhysicalPrimaryEta16++;
 855         }
 856         if(eta>-3.7 && eta<-1.7) {
 857           nChargedMCEtam37tm17++;
 858           if(isPrim) nChargedMCPrimaryEtam37tm17++;
 859           if(isPhysPrim) nChargedMCPhysicalPrimaryEtam37tm17++;
 860         }
 861         if(eta> 2.8 && eta< 5.1) {
 862           nChargedMCEta28t51++;
 863           if(isPrim) nChargedMCPrimaryEta28t51++;
 864           if(isPhysPrim) nChargedMCPhysicalPrimaryEta28t51++;
 865         }
 866       }
 867     }
 868     Int_t nChargedMC=nChargedMCEta10;
 869     Int_t nChargedMCPrimary=nChargedMCPrimaryEta10;
 870     Int_t nChargedMCPhysicalPrimary=nChargedMCPhysicalPrimaryEta10;
 871
 872
 873     // Compute the Nch weights (reference is Ntracklets within |eta|<1.0)
 874     if(fUseNchWeight){
 875       Double_t tmpweight = 1.0;
 876       if(nChargedMCPhysicalPrimary<=0) tmpweight = 0.0;
 877       else{
 878         Double_t pMeas = fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nChargedMCPhysicalPrimary));
 879         //      printf(" pMeas=%2.2f  and histo MCNch %s \n",pMeas,fHistoMCNch);
 880         Double_t pMC = fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nChargedMCPhysicalPrimary));
 881         tmpweight = pMC>0 ? pMeas/pMC : 0.;
 882       }
 883       nchWeight *= tmpweight;
 884       AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,nchWeight));
 885     }
 886
 887     // Now recompute the variables in case another MC estimator is considered
 888     // Int_t nChargedMCEta10 = nChargedMC;
 889     // Int_t nChargedMCEta16 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.6,1.6);
 890     // Int_t nChargedMCEta05 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.5,0.5);
 891     // Int_t nChargedMCEta03 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-0.3,0.3);
 892     // Int_t nChargedMCEtam37tm17 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-3.7,-1.7);
 893     // Int_t nChargedMCEta28t51 = AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,2.8,5.1);
 894     // Int_t nChargedMCPrimaryEta10 = nChargedMCPrimary;
 895     // Int_t nChargedMCPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.6,1.6);
 896     // Int_t nChargedMCPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.5,0.5);
 897     // Int_t nChargedMCPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-0.3,0.3);
 898     // Int_t nChargedMCPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-3.7,-1.7);
 899     // Int_t nChargedMCPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,2.8,5.1);
 900     // Int_t nChargedMCPhysicalPrimaryEta10 = nChargedMCPhysicalPrimary;
 901     // Int_t nChargedMCPhysicalPrimaryEta16 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.6,1.6);
 902     // Int_t nChargedMCPhysicalPrimaryEta05 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.5,0.5);
 903     // Int_t nChargedMCPhysicalPrimaryEta03 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-0.3,0.3);
 904     // Int_t nChargedMCPhysicalPrimaryEtam37tm17 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-3.7,-1.7);
 905     // Int_t nChargedMCPhysicalPrimaryEta28t51 = AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,2.8,5.1);
 906     if(fMCPrimariesEstimator==kEta10to16){
 907       nChargedMC = nChargedMCEta16 - nChargedMCEta10;
 908       nChargedMCPrimary = nChargedMCPrimaryEta16 - nChargedMCPrimaryEta10;
 909       nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta16 - nChargedMCPhysicalPrimaryEta10;
 910     } else if(fMCPrimariesEstimator==kEta05){
 911       nChargedMC = nChargedMCEta05;
 912       nChargedMCPrimary = nChargedMCPrimaryEta05;
 913       nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta05;
 914     } else if(fMCPrimariesEstimator==kEta03){
 915       nChargedMC = nChargedMCEta03;
 916       nChargedMCPrimary = nChargedMCPrimaryEta03;
 917       nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta03;
 918     } else if(fMCPrimariesEstimator==kEtaVZERO){
 919       nChargedMC = nChargedMCEtam37tm17 + nChargedMCEta28t51;
 920       nChargedMCPrimary = nChargedMCPrimaryEtam37tm17 + nChargedMCPrimaryEta28t51;
 921       nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEtam37tm17 + nChargedMCPhysicalPrimaryEta28t51;
 922     } else if(fMCPrimariesEstimator==kEtaVZEROA){
 923       nChargedMC = nChargedMCEta28t51;
 924       nChargedMCPrimary = nChargedMCPrimaryEta28t51;
 925       nChargedMCPhysicalPrimary = nChargedMCPhysicalPrimaryEta28t51;
 926     }
 927
 928     // Here fill the MC correlation plots
 929     if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
 930       fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary,nchWeight);
 931     }
 932
 933     fHistNtrVsNchMC->Fill(nChargedMC,countMult,nchWeight);
 934     fHistNtrCorrVsNchMC->Fill(nChargedMC,countCorr,nchWeight);
 935
 936     fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countMult,nchWeight);
 937     fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countCorr,nchWeight);
 938
 939     fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countMult,nchWeight);
 940     fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countCorr,nchWeight);
 941
 942     fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary,nchWeight);
 943   }
 944
 945   Int_t nCand = arrayCand->GetEntriesFast();
 946   Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
 947   Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 948   Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
 949   Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 950
 951   // pdg of daughters needed for D* too
 952   UInt_t pdgDgDStartoD0pi[2]={421,211};
 953
 954   Double_t aveMult=0.;
 955   Double_t nSelCand=0.;
 956   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 957     AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
 958     AliAODRecoCascadeHF *dCascade = NULL;
 959     if(fPdgMeson==413) dCascade = (AliAODRecoCascadeHF*)d;
 960
 961     fHistNEvents->Fill(7);
 962     if(fUseBit && !d->HasSelectionBit(selbit)){
 963       fHistNEvents->Fill(8);
 964       continue;
 965     }
 966
 967     Double_t ptCand = d->Pt();
 968     Double_t rapid=d->Y(fPdgMeson);
 969     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 970     if(!isFidAcc) continue;
 971
 972     Int_t labD=-1;
 973     if(fReadMC) {
 974       if(fPdgMeson==413){
 975         labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
 976       } else {
 977         labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
 978       }
 979       FillMCMassHistos(arrayMC,labD, countMult,nchWeight);
 980     }
 981
 982     Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
 983     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 984     if(passTopolCuts==0) continue;
 985     nSelectedNoPID++;
 986     fHistNEvents->Fill(9);
 987     if(passAllCuts){
 988       nSelectedPID++;
 989       fHistNEvents->Fill(10);
 990     }
 991     Double_t multForCand = countCorr;
 992
 993     if(fSubtractTrackletsFromDau){
 994       // For the D* case, subtract only the D0 daughter tracks <=== FIXME !!
 995       AliAODRecoDecayHF2Prong* d0fromDstar = NULL;
 996       if(fPdgMeson==413) d0fromDstar = (AliAODRecoDecayHF2Prong*)dCascade->Get2Prong();
 997
 998       for(Int_t iDau=0; iDau<nDau; iDau++){
 999         AliAODTrack *t = NULL;
1000         if(fPdgMeson==413){ t = (AliAODTrack*)d0fromDstar->GetDaughter(iDau); }
1001         else{ t = (AliAODTrack*)d->GetDaughter(iDau); }
1002         if(!t) continue;
1003         if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
1004           if(multForCand>0) multForCand-=1;
1005         }
1006       }
1007     }
1008     Bool_t isPrimary=kTRUE;
1009     Double_t trueImpParXY=9999.;
1010     Double_t impparXY=d->ImpParXY()*10000.;
1011     Double_t dlen=0.1; //FIXME
1012     Double_t mass[2];
1013     if(fPdgMeson==411){
1014       mass[0]=d->InvMass(nDau,pdgDau);
1015       mass[1]=-1.;
1016       if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
1017     }else if(fPdgMeson==421){
1018       UInt_t pdgdaughtersD0[2]={211,321};//pi,K
1019       UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
1020       mass[0]=d->InvMass(2,pdgdaughtersD0);
1021       mass[1]=d->InvMass(2,pdgdaughtersD0bar);
1022       if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
1023     }else if(fPdgMeson==413){
1024       // FIXME
1025       mass[0]=dCascade->DeltaInvMass();
1026       mass[1]=-1.;
1027       if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.0015) nSelectedInMassPeak++; //1 MeV for now... FIXME
1028     }
1029     for(Int_t iHyp=0; iHyp<2; iHyp++){
1030       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
1031       Double_t invMass=mass[iHyp];
1032       Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
1033
1034       if(fReadMC){
1035
1036         if(fPdgMeson==413){
1037           labD = dCascade->MatchToMC(fPdgMeson,421,(Int_t*)pdgDgDStartoD0pi,(Int_t*)pdgDau,arrayMC);
1038         } else {
1039           labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
1040         }
1041
1042         Bool_t fillHisto=fDoImpPar;
1043         if(labD>=0){
1044           AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1045           Int_t code=partD->GetPdgCode();
1046           if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
1047           if(code<0 && iHyp==0) fillHisto=kFALSE;
1048           if(code>0 && iHyp==1) fillHisto=kFALSE;
1049           if(!isPrimary){
1050             if(fPdgMeson==411){
1051               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
1052             }else if(fPdgMeson==421){
1053               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
1054             }else if(fPdgMeson==413){
1055               trueImpParXY=0.; /// FIXME
1056             }
1057             Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
1058             if(fillHisto && passAllCuts){
1059               fHistMassPtImpPar[2]->Fill(arrayForSparse);
1060               fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
1061             }
1062           }else{
1063             if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
1064           }
1065         }else{
1066           if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
1067         }
1068
1069         if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
1070         if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
1071
1072       }
1073
1074       if(fPdgMeson==421){
1075         if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
1076         if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
1077       }
1078
1079       fPtVsMassVsMultNoPid->Fill(multForCand,invMass,ptCand);
1080
1081       if(fPdgMeson==421){
1082         if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
1083         if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
1084       }
1085       if(passAllCuts){
1086         aveMult+=multForCand;
1087         nSelCand+=1.;
1088         fPtVsMassVsMult->Fill(multForCand,invMass,ptCand,nchWeight);
1089         fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand,nchWeight);
1090         // Add separation between part antipart
1091         if(fPdgMeson==411){
1092           if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
1093           else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
1094         }else if(fPdgMeson==421){
1095           if(passAllCuts&1) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
1096           if(passAllCuts&2) fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
1097         }else if(fPdgMeson==413){
1098           if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
1099           else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
1100         }
1101
1102         if(fDoImpPar){
1103           fHistMassPtImpPar[0]->Fill(arrayForSparse);
1104         }
1105
1106       }
1107
1108     }
1109   }
1110   if(fSubtractTrackletsFromDau && nSelCand>0){
1111     aveMult/=nSelCand;
1112     fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)(aveMult+0.5001));
1113   }else{
1114     fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countCorr);
1115   }
1116
1117
1118   fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
1119   fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
1120   fHistNtrUnCorrEvSel->Fill(countMult,nchWeight);
1121   fHistNtrCorrEvSel->Fill(countCorr,nchWeight);
1122   if(nSelectedPID>0) {
1123     fHistNtrUnCorrEvWithCand->Fill(countMult,nchWeight);
1124     fHistNtrCorrEvWithCand->Fill(countCorr,nchWeight);
1125     if(fKeepCorrPlots){
1126       fHistNtrEta16vsNtrEta1EvWithCand->Fill(countTreta1,countTreta16);
1127       fHistNtrEta05vsNtrEta1EvWithCand->Fill(countTreta1,countTreta05);
1128       fHistNtrEta03vsNtrEta1EvWithCand->Fill(countTreta1,countTreta03);
1129       fHistNtrEtaV0AvsNtrEta1EvWithCand->Fill(countTreta1,vzeroMultA);
1130       fHistNtrEtaV0MvsNtrEta1EvWithCand->Fill(countTreta1,vzeroMult);
1131       fHistNtrEtaV0AvsV0AEqEvWithCand->Fill(vzeroMultA,vzeroMultAEq);
1132       fHistNtrEtaV0MvsV0MEqEvWithCand->Fill(vzeroMult,vzeroMultEq);
1133       fHistNtrCorrEta1vsNtrRawEta1EvWithCand->Fill(countTreta1,countTreta1corr);
1134       fHistMultCorrvsMultRawEvWithCand->Fill(countMult,countCorr);
1135     }
1136   }
1137   if(nSelectedInMassPeak>0) {
1138     fHistNtrUnCorrEvWithD->Fill(countMult,nchWeight);
1139     fHistNtrCorrEvWithD->Fill(countCorr,nchWeight);
1140     if(fKeepCorrPlots){
1141       fHistNtrEta16vsNtrEta1EvWithD->Fill(countTreta1,countTreta16);
1142       fHistNtrEta05vsNtrEta1EvWithD->Fill(countTreta1,countTreta05);
1143       fHistNtrEta03vsNtrEta1EvWithD->Fill(countTreta1,countTreta03);
1144       fHistNtrEtaV0AvsNtrEta1EvWithD->Fill(countTreta1,vzeroMultA);
1145       fHistNtrEtaV0MvsNtrEta1EvWithD->Fill(countTreta1,vzeroMult);
1146       fHistNtrEtaV0AvsV0AEqEvWithD->Fill(vzeroMultA,vzeroMultAEq);
1147       fHistNtrEtaV0MvsV0MEqEvWithD->Fill(vzeroMult,vzeroMultEq);
1148       fHistNtrCorrEta1vsNtrRawEta1EvWithD->Fill(countTreta1,countTreta1corr);
1149       fHistMultCorrvsMultRawEvWithD->Fill(countMult,countCorr);
1150     }
1151   }
1152
1153   PostData(1,fOutput);
1154   PostData(2,fListCuts);
1155   PostData(3,fOutputCounters);
1156
1157   return;
1158 }
1159
1160 //________________________________________________________________________
1161 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
1162   // Histos for impact paramter study
1163   // mass . pt , impact parameter , decay length , multiplicity
1164
1165   Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
1166   Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
1167   Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
1168
1169   fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
1170                                         "Mass vs. pt vs.imppar - All",
1171                                         5,nbins,xmin,xmax);
1172   fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
1173                                         "Mass vs. pt vs.imppar - promptD",
1174                                         5,nbins,xmin,xmax);
1175   fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
1176                                         "Mass vs. pt vs.imppar - DfromB",
1177                                         5,nbins,xmin,xmax);
1178   fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1179                                         "Mass vs. pt vs.true imppar -DfromB",
1180                                         5,nbins,xmin,xmax);
1181   fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
1182                                         "Mass vs. pt vs.imppar - backgr.",
1183                                         5,nbins,xmin,xmax);
1184   for(Int_t i=0; i<5;i++){
1185     fOutput->Add(fHistMassPtImpPar[i]);
1186   }
1187 }
1188
1189 //________________________________________________________________________
1190 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
1191 {
1192   // Terminate analysis
1193   //
1194   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
1195
1196   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1197   if (!fOutput) {
1198     printf("ERROR: fOutput not available\n");
1199     return;
1200   }
1201
1202   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1203   if(!fHistNEvents){
1204     printf("ERROR: fHistNEvents not available\n");
1205     return;
1206   }
1207   printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
1208
1209   return;
1210 }
1211 //_________________________________________________________________________________________________
1212 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1213   //
1214   // checking whether the mother of the particles come from a charm or a bottom quark
1215   //
1216
1217   Int_t pdgGranma = 0;
1218   Int_t mother = 0;
1219   mother = mcPartCandidate->GetMother();
1220   Int_t istep = 0;
1221   Int_t abspdgGranma =0;
1222   Bool_t isFromB=kFALSE;
1223   //  Bool_t isQuarkFound=kFALSE;
1224   while (mother >0 ){
1225     istep++;
1226     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1227     if (mcGranma){
1228       pdgGranma = mcGranma->GetPdgCode();
1229       abspdgGranma = TMath::Abs(pdgGranma);
1230       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1231         isFromB=kTRUE;
1232       }
1233       //      if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1234       mother = mcGranma->GetMother();
1235     }else{
1236       AliError("Failed casting the mother particle!");
1237       break;
1238     }
1239   }
1240
1241   if(isFromB) return 5;
1242   else return 4;
1243 }
1244
1245
1246
1247 //____________________________________________________________________________
1248 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
1249   // Get Estimator Histogram from period event->GetRunNumber();
1250   //
1251   // If you select SPD tracklets in |eta|<1 you should use type == 1
1252   //
1253
1254   Int_t runNo  = event->GetRunNumber();
1255   Int_t period = -1;   // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1256                        // pPb: 0-LHC13b, 1-LHC13c
1257   if (fisPPbData) {
1258       if (runNo>195343 && runNo<195484) period = 0;
1259       if (runNo>195528 && runNo<195678) period = 1;
1260       if (period < 0 || period > 1) return 0;
1261   }
1262    else {
1263       if(runNo>114930 && runNo<117223) period = 0;
1264       if(runNo>119158 && runNo<120830) period = 1;
1265       if(runNo>122373 && runNo<126438) period = 2;
1266       if(runNo>127711 && runNo<130841) period = 3;
1267       if(period<0 || period>3) return 0;
1268
1269 }
1270
1271   return fMultEstimatorAvg[period];
1272 }
1273
1274 //__________________________________________________________________________________________________
1275 void AliAnalysisTaskSEDvsMultiplicity::CreateMeasuredNchHisto(){
1276   // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1277   //
1278   // for Nch  > 70 the points were obtainedwith a double NBD distribution
1279   // TF1 *fit1 = new TF1("fit1","[0]*(TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1])))*(TMath::Power(([2]/[1]),x))*(TMath::Power((1+([2]/[1])),-x-[1]))"); fit1->SetParameter(0,1.);// normalization constant
1280   // fit1->SetParameter(1,1.63); // k parameter
1281   // fit1->SetParameter(2,12.8); // mean multiplicity
1282   Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1283                         10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1284                         20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1285                         30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1286                         40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1287                         50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1288                         60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1289                         80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1290                         100.50,102.50};
1291   Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1292                     0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1293                     0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1294                     0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1295                     0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1296                     0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1297                     0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1298                     0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1299                     0.00000258};
1300
1301   if(fHistoMeasNch) delete fHistoMeasNch;
1302   fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1303   for(Int_t i=0; i<81; i++){
1304     fHistoMeasNch->SetBinContent(i+1,pch[i]);
1305     fHistoMeasNch->SetBinError(i+1,0.);
1306   }
1307 }
1308
1309 //__________________________________________________________________________________________________
1310 void AliAnalysisTaskSEDvsMultiplicity::FillMCMassHistos(TClonesArray *arrayMC, Int_t labD, Int_t countMult,Double_t nchWeight)
1311 {
1312   //
1313   // Function to fill the true MC signal
1314   //
1315
1316   if(labD>=0){
1317     AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
1318     Double_t mass = partD->M();
1319     Double_t pt = partD->Pt();
1320     fPtVsMassVsMultMC->Fill(countMult,mass,pt,nchWeight);
1321   }
1322
1323 }