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