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