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