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