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