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