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