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